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VOLTERRA  TRANSFER  FUNCTIONS  FROM  PULSE  TESTS  FOR  MILDLY  NONLINEAR  CHANNELS 


SUMMARY 


The  objective  of  this  effort  is  the  development  of  the  analytical  and 
experimental  procedures  necessary  for  equipment  EMC  characterization  in  terms 
of  nonlinear  (Volterra)  transfer  functions.  The  primary  emphasis  is  on  the 
development  and  exploitation  of  efficient  network  identification  methods  for 
determining  the  nonlinear  circuit  transfer  function  models  of  the  RF  emitter 
ports,  the  RF  coupling  paths,  and  the  RF  receptor  ports  of  a  system  level 
model  used  for  intrasystem  interference  analysis  and  prediction.  A  secondary 
emphasis  is  on  the  specification,  synthesis,  and  design  of  suitable  inter¬ 
ference  compensation  transfer  function  which  when  appropriately  combined  with 
the  identified  nonlinear  port  transfer  functions,  can  effectively  suppress 
selective,  undesired  nonlinear  responses  to  acceptable  levels,  consistent 
with  some  overall  system  performance  measure.  This  report  describes  a  method 
for  determining  linear  and  quadrative  Volterra  transfer  functions  from  pulse 
tests.  The  method  involves  two  transient  tests  in  the  laboratory,  followed 
by  data  analysis  by  the  computer,  and  is  readily  extendable  to  cubic  and 
even  higher  order  transfer  functions  where  the  computations  grow  very  rapidly. 

I .  INTRODUCTION 

One  of  the  main  problems  encountered  in  the  design  and  testing  of 
frequency-division-multiplexed  systems,  e.g.,  analog  coaxial  cable,  is  the 
determination  of  intermodulation  distortion  arising  from  the  repeater  ampli¬ 
fiers.  In  these  multichannel  systems,  many  intermodulation  products  generated 
in  an  amplifier  often  have  the  same  frequency  and  result  in  particularly  high 
interference  levels  at  certain  frequencies.  Furthermore,  when  the  number  of 
amplifiers  used  in  the  system  is  large,  as  is  usually  the  case,  certain 
distortion  products  (or  nonlinear  mixes)  are  generated  and  propagate  in  phase 
along  the  line  with  little  or  no  losses.  In  addition  other  distortion  products 
may  also  arise  in  the  ancillary  networks  used  in  conjunction  with  the  repeater 
amplifiers;  for  example,  directional  filters,  power  filters,  compandors,  and 
couplers  may  combine  and  generate  intermodulation  products.  These 
products  may  increase  to  an  exceedingly  high  level  due  to  the  accumulation 
of  these  distortion  waveforms  along  the  line.  Therefore  it  is  important 
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in  testing  such  systems  to  make  certain  that  the  distortion  lies  within 


acceptable  Limits.  The  Volterra  series  expansion  [ 1 ) — [ 3 ]  permits  description 
of  the  nonlinear  system  in  a  compact  form  and,  in  turn,  enables  expression 
of  the  distortion  in  terms  of  the  multivariable  transfer-functions  [4]-[6J. 

In  this  paper  we  describe  a  method  for  determining  the  Volterra  transfer- 
functions  H^(Sj)  and  H^Cs^.s^)  from  pulse  tests.  The  method  involves  two 
transient  tests  in  the  laboratory,  followed  by  analysis  of  the  computer.  The 
latter  consists  of  (a)  pole  determination  using  the  pencil-of-functions 
method  [7]— [9J,  and  (b)  computation  of  the  residues  by  a  least-squares 
technique.  Advantages  of  the  method  include  the  rapidity  of  the  laboratory 
tests,  as  contrasted  with  traditional  frequency-scan  approaches,  and  the 
explicit  determination  of  the  transfer  functions.  Furthermore,  the  method 
is  readily  extendible  to  H^(Sj,s,,s^)  and  even  to  higher  order  transfer 
functions,  although  the  computations  grow  very  rapidly  for  these  cases. 

The  structure  of  the  paper  is  as  follows.  In  section  II  we  give  a 
brief  description  of  the  pencil-of-functions  method  for  the  linear  transfer 
function  case.  Section  III  discusses  the  nature  of  the  Volterra  response 
of  a  2-variable  nonlinear  system.  The  identification  of  such  nonlinear 
systems  is  discussed  in  Section  IV  and  some  examples  are  given  in  Section  V. 
User  guides  to  associated  computer  programs  are  given  in  Sections  VI  and  VII. 

II.  POLE-ZERO  MODELING  BY  PENCIL-OF-FUNCTIONS  METHOD 
Recorded  input,  output  responses  of  a  network  can  be  integrated,  or 
first-order  filtered,  to  yield  a  family  of  signals,  called  information  signals. 
Application  of  the  pencil-of-functions  theorem  [9]  to  this  family  yields,  in  a 
closed  form,  the  identified  parameters  1  the  network  function.  The  procedure 
for  this  pole-zero  modeling  netln.  .  is  described  below,  and  the  program  listing 
( I GRAM)  for  the  identification  routine  is  given  in  Appendix  A. 
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Linear  Identif Icatlon  Problem 

Given  the  input-output  observations 

{u(k)>,  (y(k)},  fe  -  0,  1,  ....  K-l  (3) 

arising  from  a  physical  system  (see  Fig.  1)  believed  to  be  linear,  and  of 
finite  order,  it  is  desired  to  find  a  system  model 


H(z) 


b , z  1  +  ...  +  b  z-0 
_ 1 _ n _ 

1+  a.z  ^  +  ...  +  a  z  n 
l  n 


(4) 


R.z 

i 


n 

Z 

i=l  1  -  z.z 
x 


-1 


-1 


(5) 


which  best  fits  the  observations  (in  some  sense).  A  solution  can  be  obtained 
by  use  of  the  pencil-of-f unctions  method  [7],  [8].  The  final  formula  is 


z“A[  E/D  (1-qz'1)  ]  /D 

i=l  n  1  1 


n-i 


H(z)  - 


[  E  /D  (1-qz'1)  ]  /D 

i=l 


n-i 


(6) 


where  q  is  the  pole  of  each  of  the  first-order  filters  in  a  filter  cascade  em¬ 
ployed  for  generating  a  set  of  information  signals  y_  ,  ...,  y  ,  u  ,  ...,  u 

U  n  u  n 

(see  Fig.  2).  The  numbers  D^  are  the  diagonal  cofactors  of  the  covariance 


matrix  of  the  information  signals  (omitting  u^)  and  D  =  +  ...  +  /D2f^ 

A  tutorial  example  of  the  application  of  Equation  (6)  is  given  in  Appendix  B. 
Equivalent  s  and  z  Domain  Functions 

Conversion  between  s  and  z  domain  functions  can  be  carried  out  on  the 
basis  of  impulse-invariant  or  step-invariant  transformations  [10]: 


n  K.  Imp.  Inv.  n  K. 

Z  — —  +->-  Z  


i=l  S"Si 


i=l(l-zi  z  X) 


(7a) 


Step  Inv.  n 
«-  Z 
i=l 


-1 


K 


1 

(1-Zj^  z  ) 


(7b) 
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Fig.  1  Response-pair  from 
system  under  test 


Q(z)  =  i/n-qz"1) 


Fig.  2  Generation  of  information  signals 
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where  zi  =  e  =  K^A  and  1C  =  K^( 1-z^) /s^ : A  is  the  sampling  interval.  W 
shall  use  the  step-invariant  transformation  in  this  paper.  The  program  list  in 
for  s  to  z  (STOZ)  is  given  in  Appendix  C. 

III.  VOLTERRA  RESPONSE  TO  SQUARE- PULSE  INPUTS 
A.  Volterra  Series  Representation  of  Nonlinear  Systems 

In  the  analysis  of  wide-band  amplifiers,  it  is  often  assumed  that  the 
output  signal  depends  only  on  the  input  signal  at  the  same  instant  of  time. 
The  input/output  relation  can  thus  be  expressed  with  a  power  series  expan¬ 
sion  as  follows: 


y(t) 


ajX(t) 


+  a^x 


(t) 


+  a^x 


‘(C) 


(8) 


where  x(t)  and  y(t)  denote  the  input  and  output  signals,  respectively,  and 
the  coefficients  a^  are  time-independent  constants.  In  general,  however, 
the  output  y(t)  is  also  dependent  on  the  past  input  signal.  A  generalization 
of  (8)  in  this  case  is  a  series  of  convolution  integrals 

y(t)  =  yL(t)  +  y2(t)  +  y^t)  +  ...  (9a) 

where 


V‘> 


WV 


•  *  * ,T  )x(t 

n 


Tj)x(t  -  T  2> 


•■♦x(t  -  T  )  dT,dT„***dT  (9b) 

n  1  Z  n 

and  h^  is  a  real-valued  symmetric  function  of  n  real  variables.  Expression 

(9)  is  usually  referred  to  as  the  Volterra  series.  This  representation 

shows  that  a  nonlinear  system  may  be  regarded  as  the  combination  of  a  linear 

and  a  number  of  higher  order  nonlinear  subsystems.  Each  of  these  subsystems 

is  characterized  bv  an  impulse  response  h  (T1,T_,,*,T  ),  which  is  also  called 

n  i  Z  n 

a  Volterra  kernel.  For  a  physically  realizable  system,  h^  will  have  the 
value  zero  for  all  n  whenever  any  of  its  argument  is  negative.  Also,  these 
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kernels  are  absolutely  summable  for  stability.  The  nth  order  transfer 


function  is  defined  as  the  n-fold  Laplace  transform  [11]  of  h  ,  i.e., 


H  (s.  ... ,s  )  = 
n  1  11 


f 


-(s  T ST  ) 

h  (i  ,T  ,, .1  )e  n  dT  dT  (10) 

n  1  2  n  In 


In  particular,  we  shall  call  H^(s^)  the  linear  transfer  function 
and  the  quadratic  transfer  function,  and  the  corresponding 

contributions  to  the  response  as  the  linear  and  quadratic  responses, 
respect ively . 

B .  Two- Var iab le  Volterra  Svstem 

The  two  variable  Volterra  system  is  charact erized  by  eouation  (8) 
when  it  is  terminated  at  n  =  2 ,  i.e., 
y(t)  =  y^t)  +  v  ^  ( t ) 

In  the  rest  of  the  paper  we  shall  assume  that  the  quadratic  subsystem 

H^s^.s,,)  has  the  form  shown  in  Fig.  3. 

In  testing  a  two-variable  Volterra  system,  the  following  remark  is 

very  useful.  It  is  possible  to  separate  the  linear  and  quadratic  response 

by  performing  two  measurements.  Specifically,  let  the  response  of  the 

system  to  the  inputs  x+(t)  =  x(t)  and  x  (t)  =  -x(t)  be  denoted  as  y+(t)  and 
y~(t)  respectively.  From  (9)  it  follows  that  y+(t)  =  y^Ct)  +  y ^  ( t )  and  y“(t) 
-y^t)  +  y2(t)  ,  so  that 

y  1  (t; )  =  -|[y+(t)  -  y  ( t)  ]  (11a) 

y2(t)  =  -|[y+(t)  +  y  (t)l  (lib) 

Clearly,  H^(s^)  can  Identified  from  the  pair  x(t)  and  y^(t)  by  using 
the  technique  of  Section  II.  Examples  of  successful  identification  of 
linear  kernels  are  given  in  (7]  and  [12];  the  former  includes  the  case  of  a 
wideband  RF  amplifier.  Fig.  A  is  reproduced  from  reference  [7]. 


We  shall  therefore  concentrate  on  the  identification  of  H_(s,,s0)  from 


V 


v 


H^M'  *2* 


Ha(s)  •  X 
i=1 


Ai 


(s  +  ai ) 


Hc(s)  -  X 
k=l 


(5  +  Ck) 


Fig.  3  2-Variable  system:  Linear  subsys 
and  quadratic  subsystem  H2 ( s 
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the  pair  x(t)  and  y9(t). 

C.  Souare-Pulse  Response  of  Quadratic  TF 


From  the  block-diagram  of  Fig.  3,  it  can  be  readily  shown  that 

H2(S1,S2)  =  Ha(sl)Ha(s2)Hc(sl+S2)  (12) 

so  that  the  associated  two-dimensional  response  is^ 


Y^2)(si»s2)  =  Ha(s1)Ha(s2)Hc(s1+s2)X(s1)X(s2) 


(13a) 


For  a  square  pulse  input  p(t)  =  u(t)  -  u(t-T),  the  associated  response  for 

2 

the  specific  transfer  functions  of  Fig.  3  becomes 


-s  T  -s  T 
n  A.A.Ck(l-e  1  )  (1-e  ) 


(13b) 


From  this  the  response  Y2(s)  can  be  found  at  once  by  use  of  the  theorem  of 

2 

Appendix  D.  The  inverse  transform  yields 

y2(t)  = 

n  AiAiSc  k  ik  -ait  ik  -ait  ilk  _cict 

Z  1  [d  -die  -d,  e  J  +d::i  e  1  J  +  fdJ  - 

i.j.k-i  Vj  01  1  3 


+(dj  -d5  )e  K  ]u(t) 


for  0  <  t  <  T 


n  A  A.C,  .  -(a  +a  )t’ 

I  --  J  k  [  dJk(L1-l)(Lj-l)e  1  1  + 

i.j.k-i  Vj  1 


((d^k-d3;lk)Nk+d5Li+d^ikLj-d3:’k-d3k)e  k  ]  u(c’) 


where 


for  T  <  t 


(14) 


t’  =  t  -  T 

-a,  T 

L1  =  e  1 

k  "CkT 

N  =  e 


*  Parantheses  around  2  are  used  to  emphasize  that  Y(2)  is  tke  associated 
two-dimensional  response. 


Note  that  if  H  has  m  n  poles,  then  the  index  k  runs  from  1  to  m  in  (13) 
and  (14) 


f 


lk  = 

l/c. 

o 

k 

- 

1/(Vai) 

,ijk  _ 

3 

l/(ck-a.-a.) 

II 

*r-> 

aj/ck(ck-aj) 

,ijk  _ 

a,/(c, -a.  )(c 

(15) 


IV.  IDENTIFICATION  OF  USING  PULSE  INPUTS 

In  this  section  we  deal  with  the  central  problem  of  the  paper  i.e., 

identification  of  H2(si’s2^  from  the  Volterra  response  y2(t).  The  problem 

is  considered  in  two  parts.  First  the  estimation  of  poles  of  H  and  H  is 

a  c 

considered.  Next,  the  residues  are  estimated. 

A.  Identification  of  Poles 

A.l  Poles  from  Response  over  0  £  t  £  T  (Segment  1) 

For  the  symmetric  quadratic  subsystem  the  response  to  the  square-pulse 
p(t)  =  u(t)  -  u(t-T)  was  shown  to  be  given  by  (14).  Over  the  time  interval 
0  <  t  <  T  this  response  can  be  visualized  as  the  unit-step  (at  t=0)  response  of  a 
linear  system  shown  in  Fig.  5(a)  where  A  a^  +  a^  and  the  residues  P^, 
and  R^  are  defined  according  to  the  first  part  of  (14). 

Clearly,  the  pencil-of- f unct ions  method  [7], [9]  can  be  used  to  determine 
a.  i  =  1,  . .  . ,  n 


a.  . 
ij 


i=l, . . . ,n; 
k  1,...., m 


j=l, • • ■ ,i 


Note  that  in  general  the  total  number  of  poles  is 
n(n+!) 

N  =  n  +  — ^ - +  m 


(16) 
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If  the  quadratic  is  completely  symmetric,  i.e.,  if  H  =  If  .  then  the 
-  ■  '  a  c 

total  number  of  poles  is  again  given  by  (16);  however,  the  poles  a^  occur 
with  a  multiplicity  of  2. 

A. 2  Poles  from  Response  over  T  £  t  (Segment  2) 

Over  the  time  interval  T  £  t  the  quadratic  response  is  given  by  the 
second  part  of  (14),  and  can  be  visualized  as  the  unit  impulse  (at  t=T)  response 
of  a  linear  system  as  shown  in  Fig.  5(b) 

Again  the  pencil-of-functions  method  can  be  applied  to  determine 


cfc  k  =  1, . . .  ,m 

Note  that  now  the  total  number  of  poles  is 

_  n (n+l) 

N  =  — -  +  m 


(17) 


A.  3  Remarks 

1.  The  test  engineer  has  a  choice  here  between  the  use  of  the  two 

segments  of  At  this  time  it  appears  that  the  use  of  Segment  2  is 

preferable,  since  the  dimensionality  of  identification  is  lower  in  this 
case  (as  evidenced  from  a  comparison  of  (17)  with  (16)),  although  quite 
extensive  experimentation  is  necessary  to  make  a  definitive  statement.  For 
brevity  of  the  paper  we  will  assume  that  Segment  2  is  used  for  identification, 
and  will  give  details  only  for  this  case. 

2.  Since  the  identified  values  a  and  c  will  generally  not  coincide  with  tin 

1  j  K 

true  values,  it  is  necessary  to  isolate  the  poles  by  visual  inspection  or 
by  a  suitable  computer  routine.  For  example,  if  these  numbers 
for  the  case  (n,m)  =  (2,1)  are 

2.1,  4.15,  9.8,  5.95 


then 
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A 

311  " 

2.1 

a12 

4.15 

a22 

5.95 

A 

C1  " 

9.8 

and 

we  may 

take 

ai  " 

1.05 

A 

a2  ■ 

3.1 

h  = 

9.8 

B. 

Identification  of  Residues 

After 

the  poles  have  been  determined. 

n  i  m 

y2(t) 

(18a) 


where  are  defined  in  accordance  with  (14)  (and  are  now  known  since  they 

are  functions  of  the  poles  and  the  pulse  width  T) .  Further,  by  defining 

(18b) 


ev  "  AiAjCk 


8V(U) 


fijk(uA) 


v»i  +  j+  k-  2 
y  =  0, . . . ,M-1 


(18c) 


obtain  the  following  set  of  simultaneous  equations 


Y  =  G  E  (19) 

where 

I  =  [y 2 (0)  y2(A) . y2((M-l)A)]T 

G  -  lgv(p)]  M  x  I  matrix 


Note  that  M  denotes  the  number  of  time  samples  used  in  (18)  and  (19), 

and  I  denotes  the  number  of  unknown  residue-products.  If  M  is  chosen  equal 

3 

If  some  of  the  poles  are  complex  (occuring,  in  conjugate  pairs),  the 
associated  residues  are  also  complex.  In  such  cases,  it  is  possible 
to  equate  real  (and  imaginary)  parts  on  both  sides  of  the  equation  (19)  to 
obtain  real  coefficient  equations 
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i 


f 


f 


to  I  then,  we  obtain  the  solution 


E  =  G'1  V  (20) 

More  generally,  if  M  >  I,  then  we  obtain  the  leas t-squares  solution  [13] 

I  =  (GTG)'L  Gty  (21) 

Finally,  the  equations  (18b)  can  be  solved  straightforwardly  to  yield  the 
residues  and  C^.  Because  of  the  homogeneity  of  the  relations  any  one  of 
tne  real  residues  (or  the  magnitude  of  a  complex  residue)  can  be  taken  as  1. 

V.  EXAMPLES 

Two  computer-generated  examples  will  be  presented.  The  first  is  a 
simple  case  with  n=m=l  and  is  intended  to  clearly  present  the  details  of 
the  procedure.  The  second  is  a  more  complicated  case,  representing  a  some¬ 
what  realistic  channel.  It  has  a  linear  transfer  function  characterized 
by  a  6th  order  butterworth  filter,  and  a  quadratic  transfer  function  with 
n=3  and  m=l. 

Example  1 

Consider  that  the  response  y  (t)  of  a  quadratic  subsystem,  with  H  (s)= 

2  a 

l/(s+l)  and  R  (s)=l/ (s+0. 5) ,  to  a  square  pulse  of  one  second  duration  has 
c 

been  measured. 

The  second  segment  of  the  response  (y2(t),  1  <  t  <  2)  is  sampled  with 
A=0.02  sec.  and  used  for  identification.  To  first  determine  the  parameters 

of  the  equivalent  linear  system,  a  unit  sample  pulse  input  is  assumed  and 
<j=0.8  is  taken  for  the  seneration  of  the  information  signals.  The  Cram 
matrix  of  the  correlation  signals  turns  out  to  be 


>.  333030+01 
i.  4  10*30+02 
1334  30+03 

0. 19507D+03 

0. 9032SD+03 

0.423710+04 

let  /7D  +  01 

-  O.  4&jC4C'+01 

-O. 123670+02 

0. 277730+01 

— 

'.91  3  3  30+01 

-O.  j 3616li ♦02 

-0. 1 16230+03 

0.771600+01  0.331510+02 

The  normalized  cofactor  square-roots  (/D\//n^)  are 

1.0000  0.35931  0.032129  0.22357  0.026502 

I3v  use  of  (6)  and  (7h)  the  negatives  of  the  poles  are  computed  to  be^ 


The  reason  to  associate  1.9997  with  aj^  is  that  when  segment  1  is  analyzed, 
a  pole  at  0.99°6  also  shows  up. 
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1.9997 


all  = 
c  =  0.5002 

so  that 

al  =  0.9998 
Cl  =  0.5002 

Note  that  the  waveform-fit  error  in  modeling  the  poles  of  the  equivalent  linear 
system  was  almost  negligible  (normalized  mean-square  error  ■  10~7)  pointing  to 
the  success  of  pulse  testing  approach. 

Now,  using  these  poles  we  find 


y2(t’)  =  A1A1C1  (-0.2664375  e"1*  "97t '+0. 4119195e“°‘ 5002t  ’  ) 

Using  a  single  point  for  numerator  computation,  y2(t’=0)  =  0.292924,  we  obtain 
A1A1C1  *  2.013,  so  that5 


1 

2.013 


Example  2 


Linear  TF;  Sixth  order  lowpass  buttervorth  filter  with  cutoff  f  =  10  MHz 

c 


6. 1528908(10^) 


Hl(s)  “  ,6 


[su  +  2.4276(108)s5  +  2.9467(1016)s4  +  2.2676(1024)s3  + 


T?  ?  19  46 

1.1633(10  )s  +  3. 78358(10 +  6.1528908(10  )] 


6.1528908(1010) 


[A6  +  242. 76A5  +  294671*  +  2267581A3  +  1. 6133(10®) A2 


q  in 

+  3.78358(10  )A  +  6. 1528908( 10  ) ] 


'a  better  value  A^A^C^  =  2.0007  is  obtained  using  50  points  and  formula  (21). 
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Where  X  =  (10  8)s  (units  of  Mrad./s). 


so.sao13) 


Quadratic  TF  H  (s)=H  (s)=  — - - - - — 

C  sZ  +  (10  )s  +  25.25  (1014) 


X  +  10X  +  2525 

The  above  2-variable  system  was  excited  by  a  square  pulse  p(t)  of 
duration  T  =  1  sec  and  the  response  y+(t)  recorded.  The  system  was  next  excited 
by  the  opposite  polarity  pulse  -p(t)  and  the  corresponding  response  y  (t) 
also  recorded.  From  these  responses  we  obtained  the  linear  response  y^(t) 
and  the  quadratic  response  y2(t)  by  use  of  (11). 

The  results  of  identification  are  given  below: 

Estimated  Linear  Transfer  Function  - 


15.2  X  +  6.1523(10  u) 

[X6  +  242. 76X5  +  29467X4  +  2267593X3  +  1. 6133 (108)X2  +  3.784(109)X 


+  6.15237(10  U)  ] 


Estimated  Quadratic  Transfer  Function 


-  -  0.003  X  +  505.07 

H  (X)  =  H  (X)=  — - 

X  +9.9998X  +  2524.8 


VI.  USER  GUIDE  FOR  THE  COMPUTER  PROGRAM  '  IGRAM' 

The  computer  program  IGRAM  has  two  functional  parts. 
The  first  deals  with  the  reading,  or  internal  generation, 
of  the  response  pair  of  the  system.  The  second  pertains 
to  identification  of  the  system.  This  second  part  also 
contains  the  noise  correction  facility  which  is  useful 
when  the  data  are  corrupted  by  noise.  On  an  optional 
basis  the  identified  z-domain  model  can  also  be  converted 
to  the  s-domain.  Output  from  the  program  is  in  two  forms 
a)  the  printed  output  consisting  of  the  identified  model 
and  the  error  of  fit,  and  b)  certain  plot  files  which  may 
be  used  for  plotting  purposes. 

The  program  uses  the  following  subroutines 

BUILDA  Constructs  the  transformation  matrix  A 
to  help  implement  (6).  (see  reference  [7]) 
BUILDZ  Constructs  the  covariance  matrix  of  the 
noise  signals  (assuming  unit  variance  at 
the  input  of  the  processing  filters) 
CORUPT  Adds  noise  to  the  system  response  if  the 
variable  'VAR'  is  non-zero.  Note  that 
this  is  normally  used  in  simulation  mode, 
i.e.,  when  a  known  system  function  H(z) 
is  used  and  a  response  pair  is  generated 
within  the  program. 

ERROR  Calculates  the  percent  mean  square  error 
and  percent  RMS  error;  the  latter  is  de¬ 
fined  as 


TOO 
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FILLV 


Generates  the  input  waveform  according  to 
the  input  option  selected.  These  waveforms 
range  from  a  simple  step  to  a  rectangular 
pulse  .  and  from  a  continuous  sine  wave 
to  a  sinusoidal  burst-  Also  included  are 
certain  special  waveforms  such  as  a  chirp 
and  triangular  pulse.  (See  page  22  for 
detail s. ) 


FIX  Performs  noise  correction  on  the  Gram 

matrix  G  to  enable  estimation  of  a  more 

reliable  system  model.  Estimated  noise 

variance  is  * 

2  1 
CT  s 

- =!— - 

G  ®W 

where  W  is  the  covariance  matrix  of 
noise  vector.  The  estimated  gram  matrix 

1S  2 

F  =  G  -  <t  W 


GKRDCT 


GRAMI 


IZTOS 


POLCON 


This  routine  finds  the  cofactors  and/or 
the  inverse  of  a  square  matrix.  It  also 
calculates  the  the  denominator  param¬ 
eters  through  pencil-of-f unctions  method. 
Performs  the  Pencil-Of-Functions  technique 
for  identification.  It  calls  several 
subroutines,  notably  GKRDCT  and  BUILDA. 
Separates  the  numerator  and  denominator 
parameters  and  collects  them  into  two 
vectors.  Also  it  calls  ZTOS  for  z 
domain  to  s  domain  conversion 
Constructs  the  polynomial  corresponding 
to  a  known  set  of  roots. 


*The  matrix  DOT  product  A9B  means  l  z 
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POLRT 


Computes  the  real  and  complex  roots  of  a 
real  polynomial.  Limited  to  a  36th 
(or  lower)  order  polynomial. 

PRCVEC  Prints  a  double  precision  complex  vector 
in  the  form  of  a+jb,  where  a  is  the  real 
part  and  b  is  the  imaginary  part. 

PR MAT  Prints  a  double  precision  two  dimensional 

array. 

PRVEC  Prints  a  double  precision  one  dimensional 
array. 

RESPON  determines  the  response  of  a  discrete 

time  transfer  function  H(z)  to  an  input 
sequence  v(k).The  coefficients  of  H(z)  are 
in  the  vector  GAMMA. 

ZTOS  Converts  the  z-domain  model  H(z)  to  an 

equivalent  s-domain  model  H(s). 

INPUT  CARDS 

Card  1  Title  card 

Card  2  Subtitle  card  (containing  list  of  variables  on 
the  next  card) 

Card  3  N=  System  order  (for  simulation  purposes) 

(1215) 

NPT=  Number  of  points  (generated  by  simulation 
or  read) 

IPLT=  0  no  plots 

=  1  plot  on  printer 

=  2  fill-in  plot  files 

IS IM  Selects  type  of  mode  desired 
=  0  lab  data  mode 

=  1  simulation  mode;  i.e.,  response- 

pair  is  generated  internally 
based  on  a  given  H(z) 
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mixed  mode;  output  data  is  lab  data, 
input  is  generated  in  the  program 
INPT  Selects  type  of  input  desired  for  simula¬ 
tion  (see  page  22  for  further  details) 
NPUL,  NCY  Parameters  of  specific  inputs,  used 
only  if  ISIM  equals  0  or  2  (see  page  22 
for  further  explanation) 


If  IS IM  =  0  use  the  following  cards 
Card  4 
to 

Card  4+NPT/7  NPT  output  data  points  (7F10.0) 

Next  NPT/7  Cards  NPT  input  data  points  (7F10.0) 

If  IS IM  =  1  use  these  cards 

Card  4  Denominator  of  H(z)  transfer  function 


Card  5  Numerator  of  H(z)  transfer  function 

Both  cards  are  (7F10.0) 

If  IS IM  =2  use  this  arrangement 
Card  4 
to 

Card  4+NPT/7  NPT  output  data  points  (7F10.0) 


Next  card 
Next  card 


Next  card 
(6  15 , F10 . 0 , 
6F5.0) 


Subtitle  card 

Subtitle  card  (containing  name  of  variables 
on  the  next  card) 

N=  model  order 

NPT  =  Number  of  points  (generated  by  sim¬ 
ulation  or  read) 

ISKIP  =  Skip  interval  in  printer  plots. 

For  example,  if  ISKIP=5  every 
5th  point  of  the  data  array  is 
plotted  (thus  reducing  the  size 
of  the  plot) 
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IREM  and  IDLY  adjust  for  the  type 
of  numerator  desired  in  the  model, 
i.  e. , 


H  (z )  =  _ 

IREM  =  0 
=  1 

=  m 

IDLY  =  0 

=  1 

=  m 

An  example 


-1 

b  +  bz  + . 

0  1 

=T 

1  +  a  z  + . 

1 

no  effect 

b  =  0  is  imposed 
n 

b  =  . . .  b  * 


-n 

+  b  z 
n 

=rr 


+  a  z 
n 


0 


n-m+1  m 

no  effect 

numerator  is  multiplied  by 
numerator  is  multiplied  by 


-1 


z 


-m 

z 


If  the  model  is  desired  to  be  of  the 
form 


H  (z)  = 


b  z 
1 


-2 


1 


- =1 -  -2 

+  a  z  +  a  z 
1  2 


then  N=2  and  IREM  and  IDLY  should  be 
set  to  1. 


DELTA 
IB  IAS 
Q 

VAR 


sampling  interval 
1  enables  bias  extraction  option 
value  of  the  processing  filter 
pole 

variance  of  the  noise  added  to 
the  simulated  response 
A  threshold  used  in  the  noise 
correction  subroutine  FIX 


DFAC 


Next  card 
(1215) 


IPR 


level  of  output  printing 
=  0  minimal  printing 
=  1  maximum  printing 


IFIX 


IPL1 


noise  correction  option 
=0  no  correction 

=1  do  noise  correction 

and  IPL2  plot  options 
=0  no  plot  files 
=  1  create  plot  data  files 


Further  explanation  of  INPT.  NTH,  NCY.  NPUL 
INPT 

Specifies  desired  input  waveform 
=  0  zero  input 
=  1  impulse  located  at  k  =  1 
=  2  step  function  (amplitude  =  0.1  *  NPUL) 

=  3  square  pulse  (width  =  NPUL) 

=  4  doublet  (total  width  =  NPUL) 

=  5  positive  triangular  pulse  (width  =  NPUL) 

=  6  triangular  pulse,  positive  and  negative 
(width  =  NPUL) 

=  7  square  wave  (period  =  NPUL) 

=  8  square  wave  burst,  followed  by  an  exponen¬ 

tially  decaying  tail  (width  =  NPUL, 
f requency=NCY/NPUL,  time  constant  =  NPT/10) 
=  9  exponential  decay  (time  constant  =  NPUL) 

=  10  impulse  train  (period  =  NPUL) 

=  11  exponentially  decaying  sine  oscillation 
(time  constant  =  NPT/5,  period  =  NPUL) 

=  12  exponentially  decaying  cosine  oscillation 
(time  constant  =  NPT/5,  period  =  NPUL) 

=13  random  signal,  sigma  =0.1  *  NPUL 
(variance  =  sigma  **2) 
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=  14  cosine  burst  (frequency  =  NCY/NPUL,  width  = 
NPUL) 

*  15  sine  burst  (frequency  =  NCY/NPUL,  width  = 
NPOL) 

=  16  cosine  chirp  (initial  frequency  =  NCY/NPUL, 
width  =  NPUL ) 

=  17  sine  chirp  (initial  frequency  *  NCY/NPUL, 
width  »  NPUL) 

NTH 

This  causes  the  program  to  subsample  the  waveforms. 
That  is,  every  nth  point  of  the  data  (measured  or 
simulated  data)  is  used,  the  rest  are  discarded. 

If  NTH=3 ,  then  every  third  point  of  the  input 
and  output  data  is  used  in  the  identification 
procedure. 

This  option  is  sometimes  useful  in  nonlinear 
systems  studies.  More  will  be  said  about  it 
elsewhere. 

NCY 

When  an  input  waveform  is  a  burst  NCY  defines 
the  number  of  oscillations  desired  in  the  time 
length  NPUL.  When  an  input  is  a  chirp  it  defines 
the  initial  number  of  cycles  per  time  duration 
of  NPUL. 

NPUL 

Any  input  that  has  a  specific  time  duration  (burst, 
chirp,  pulse,  impulse,  square  wave)  will  have  NPUL 
equal  to  the  duration  of  that  waveform. 
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VII.  USER  GUIDE  FOR  THE  COMPUTER  PROGRAM  'STOZ' 


STOZ  is  a  s-domain  to  z-domain  conversion  program. 
Given  an  s-domain  transfer  function  H(s),  it  finds 
an  equivalent  z-domain  transfer  function  H(z)  by  one 
of  several  methods,  e.g.  impulse-invariant  or  step 
invariant  transformation.  Thus,  given  the  continuous¬ 
time  description  of  a  linear  system,  it  computes  the 
equivalent  discrete- time  description. 

The  input  arrays  A  and  B  are  filled  according  to  the 
differential  equation 

a  (0 )  *y  +  a(l)*D(l,y>  +  ....  +  a (n) *D (n,y ) 

=  b  (0 ) *u  +  b  (1 )  *D  (1 ,  u)  +  _  +  b(n)*D(n,u) 

(where  D(m,f)  =  the  mth  time  Derivative  of  function 
f(t)),  or  the  transfer  function 


n 

b  +  bs  + .  +  bs 

0  1  n 

H  (  S )  =  _ 

n 

a  +  as  +  .  +  as 

0  1  n 


a  =1  (mandatory) 
n 

STOZ  returns  arrays  A  and  B  containing  the  equivalent 
discrete-time  description  stored  according  to  the 
difference  equation 

y(k)  +  a(l)*y(k-l)  +  ....  +  a(n)*y(k-n) 

=b (0 ) *u (k )  -  b ( 1 ) *u(k-l)  - _ -  b (n) *u (k-n) 


or.  the  transfer  function 


H  ( z )  = 


-1 

b  +  b  z  + 

0  1 

- _r_ 

1  +  a  z  + 

1 


+  b  z 
n 


-n 


+  a  z 
n 


-n 
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The  poles  of  the  continuous  domain  must  be  distinct  and 
non-zero  for  the  transformation  to  be  valid. 


The  program  uses  the  following  subroutines 

PCSTZ  Constructs  the  polynomial  corresponding 
to  a  known  set  of  roots;  the  kth  root  is 
omitted  if  k  is  not  equal  to  zero. 

POLRT  Computes  the  real  and  complex  roots  of  a 
real  polynomial.  Limited  to  a  36th 
(or  lower)  order  polynomial. 

PRCVEC  Prints  a  double  precision  complex  vector 
in  the  form  of  a+jb,  where  a  is  the  real 
part  and  b  is  the  imaginary  part. 

PRVEC  Prints  a  double  precision  one  dimensional 
array. 

Note  that  POLRT.  PRCVEC  and  PRVEC  are  the  same 
as  in  IGRAM,  therefore  their  listing  will  be 
omitted  in  Appendix  B. 

DATA  CARD  SET  PREPARATION 

Card  1  Title  card 

Card  2  Subtitle  card  (containing  list  of  variables 
on  the  next  card) 

Card  3  Contains  the  following  variables 
(3 15 , F10 . 0 , 15  ) 

N  =  order  of  system 

N  (maximum)  =  one  less  than  the  dimension  subscript 

IMTHD  =  0  for  the  impulse  invariant  conversion 
=  1  for  the  pulse  invariant  conversion 
=  2  for  the  trapezoidal  invariant  conversion 
=  3  for  the  logrithmic  transform  conversion 
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IPOLZ 


=  1  poles  and  zeroes  are  given  as  complex 
numbers  (real  and  imaginary);  negatives 
of  poles  and  zeros  are  read;  ie, 
for  a  pole  of  (s+2),  input  +2.0  +0.0 

(2F10.0  per  pole,  4  poles  per  card) 

=  0  denominator  and  numerator  are  read 
in  polynomial  form. 

Coefficients  ordered  from  low  to  high  degree, 
denominator  information  always  read  first. 
Highest  order  denominator  coeff  must  be  1.0 

Note:  When  pole-zero  data  is  entered  a  gain  card 
must  follow  the  last  pole  card. 

If  there  are  no  zeros,  use  blank  cards  in 
the  normal  zeros  positions. 

DELTA  =  Sampling  interval  in  seconds 
IZERO  =  1  zeros  of  z-domain  function  are  printed  out 
If  IPOLZ  =  0 

(§ff$  $)  Denominator  polynomial  coefficients 

Card  5  Numerator  polynomial  coefficients 

(8F10.0) 

If  IPOLZ  =  1 


$)  Denominator  poles  (complex  values) 

Card  5  Numerator  zeros  (complex  values) 

(8F10.0) 

Card  6  Numerator  gain  value 

(F10.0) 
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VIII  CONCLUSIONS 


Communication  channels,  both  microwave  and  satellite,  generally  exhibit 
nonlinear  characteristics.  When  frequency  division  multiplexing 
is  employed,  it  is  important  to  evaluate  the  intermodulation  distortion  using 
Volterra  models.  We  have  presented  a  modeling  procedure  which  employs  pulse 
testing  and  performs  penci 1-of-f unctions  analysis.  The  advantages  of  the 
method  are:  a)  the  use  of  a  commonly  available  signal  source;  b)  the  rapidity 
of  the  test  (only  two  transient  tests  are  required),  and  c)  the  resulting  model 

is  in  the  compact  transfer-function(s)  form.  The  success  of  the  procedure 
was  demonstrated  by  computer  generated  examples.  In  actual  application, 
however,  several  points  of  caution  must  be  exercised.  For  example,  in  the 
reverse  polarity  test,  the  amplitude  and  width  of  thepulse  must  be  closely 
duplicated.  The  development  in  thepaper  was  carried  out  under  the  assumption 
of  a  baseband  channel.  However,  with  suitable  modifications,  such  as  the 
use  of  a  carrier  frequency  burst  instead  of  a  square  pulse,  it  should  apply 
to  bandpass  channels  also.  Mathematical  details  of  this  important  extension 
remain  to  be  examined. 
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APPENDIX  A 


i 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  IGRAM  REV  10/82  V  K  JAIN 


it*********************************************************** 

• 

*  PROGRAM  NAME  'IGRAM' 

*  REVISED  10/82 


NETWORK  TRANSFER-FUNCTION  IDENTIFICATION  PROGRAM  UTILIZING 
PENCIL-OF-FU NOTIONS  METHOD 


DEVELOPED  AT  THE 
UNIVERSITY  OF  SOUTH  FLORIDA 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


FOR  ASSISTANCE 
CONTACT 
DR.  V.  K.  JAIN 
813/974-2581 


DIMENSION  X (1024 ) , V ( 1024 >, XORG ( 1024 > ,VORG(1024) ,XREC(1024) 
DIMENSION  DATA (1024 ,  2  J  , DATA2 (1024,2) , BUFF (3072) 

DIMENSION  G (20,20)  ,Z (20,20)  ,  GAMMA  (20)  ,  XLAtlDA  (20 )  ,  COEFF  (20) 
DIMENSION  GDUM (20,20) ,GEST ( 20 , 20 ) ,GGG ( 20 , 20) 

REAL*8  GDUM,GEST, GGG 
DIMENSION  H PULSE (64) 

DIMENSION  TITLE (80) 

INTEGER  LABEL < 20 ) ,CCHK 
INTEGER  ENDFIL  /  '//  '  / 

REAL*8  G,Z, GAMMA, XLAMDA, COEFF 

REAL *8  DELTA, Q, QSAV, DELSAV, AVGQ, AVGW, SUMV2 , XSAV 

EQUIVALENCE  (Z  (1,1)  ,BUFF(1)  >,(G(1,1),  BUFF  (1501)) 

EQUIVALENCE  (DATA  ( 1 , 1 )  ,  X  ( 1  > )  ,  (DATA(1,2)  ,VORG  (1) ) 

EQUIVALENCE  (XORG  (1)  ,DATA2(1,1)  )  ,  (XREC(l)  ,DATA2  (1,2)  ) 

COMMON  /GKRD/IGKR 

COMMON  /DA1/NN, IFIX 

COMMON  /DA2/IPR 

COMMON  /DA3/DFAC 

COMMON  /BIAS/IB IAS 

COMMON  /DELAY/  IDLY 

DATA  CCHK/4H  _/ 

COMMON  /WAVE/NCY 


REWIND  9 

IGKR=1 

RDEL=0 . 002 

MAX=20 

MAXPL=1024 

NPT-MAXPL 

WRITE (6,2) 
WRITE (6, 1022) 


t 
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READ (5,1021 ) TITLE 

WRITE (6,1021) (TITLE (I ) , 1-1,70) 

READ (5, 1021) TITLE 

IF ( I  PR. GE. 2 (WRITE (6, 1023) 

C  . 

c 

C  READ  LABORATORY  RECORDED  RESPONSE-PAIR, 

C  OR  GENERATE  IT  BY  SIMULATION 

c 

c 

4320  READ  (5, 1001) N,  tlPl ,  IPLT,  ISIM,  DIPT,  NPUL,  NTH , NCY 
ISIMIN-0 

IF  l  ISIM. EQ. 2) ISIMIN-1 

IF (ISIM. EQ. 2 ) ISIM=0 

I F (NCY. EQ. 0) NCY-2 

IF (NTH. EQ. 0) HTH-1 

I F  ( IS  III.  NE.  0)  GO  TO  3306 

READ  (  5, 6995)  (XORG (K) , K=1 ,MP1 ) 

C  WRITE (6,6996) (XORG <K> , K=1 ,NP1 ) 

I F  ( ISIHIN.  EQ.  0) READ  (5,6995)  (VORG  (K ) ,  K=1 , 11P1 ) 

I F  (ISIMIN.  EQ.  1)  CALL  FILLV  (VORG ,  MP1 ,  INPT,  NPUL) 

C  WRITE  (6, 3) 

C  WRITE  (6, 6996)  (VORG  (K ),  K*1 , 11P1 ) 

GO  TO  3308 
3306  CONTINUE 

MPMP2-N+N+2 

11P1-N+1 

NP2=IJ+2 

CALL  FILLV  (VORG,  HP1 ,  INPT,  NPUL) 

READ (5, 7 01) (COEFF(I) ,I=1,NP1) 

READ (5, 701) (COEFF ( I )  ,  I-NP2  ,  NPNP2 ) 

CALL  PRVECI  COEFF,  NPNP2) 

701  FORMAT (7F10.0) 

CALL  RESPON( XORG, VORG, N, COEFF, XLAMDA.NPl) 

KK-0 

DO  703  K-l  ,11P1 , NTH 
KK-KK+1 

VORG  (KK)  -VORG  (K) 

703  XORG(KK)=XORG(K) 

MP1= (MPl+NTH-1 ) /NTH 
MP2-HP1+1 

DO  705  K=MP2 , MAXPL 
VORG (K ) =0 . 0 
705  XORG (K) =0 . O 
3308  CONTINUE 
C 

4321  READ (5,1033,END=1234) (LABEL (I >, 1=1 ,20) 

IF  ( LABEL  ( 1 )  .  EQ.  ENDFIL)  GO  TO  1234 
READ (5,1021) TITLE 

READ ( 5 , 1 ) N, MP1 , ISKI P, IREM, IB I AS, IDLY, DELTA, 

+0SAV,  VAR,  DFAC 
READ (5, 1021) TITLE 
READ (5, 1001)1 PR, I FIX, IPLl , IPL2 
I F (I).  EQ.  10000) GO  TO  4320 

IF  ( IPR.  GE.  1 )  WRITE  (6,1000 )  N,  MPl ,  DELTA,  IREIl,  IB  IAS,  IDLY 
C 

IGKR=1 

IZTS-2 

IF (QSAV. EQ. 0. 0 ) QSAV=1 . 0D00 

O-OSAV 

NM1-N-1 

NP1-N+1 

NP2-N+2 

NPNPl  =  N-fN+l 

N PNP2-N+N+2 

NN-N-IREM 
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i 


I 


c 

C  ADD  NOISE  (OF  GIVEN  VARIANCE)  TO  OUTPUT,  IF  DESIRED 

C 

c 

VMAX=0.0 

xriAX=o.o 

xxsun=o.o 

23  D024 1  =  1 ,  (1P1 

I F  ( AB  S  (VORG  (I  )  )  .  GT. VMAX )  VliAX-AUS  ( VORG  ( I  )  ) 

I  F  ( ADS  ( XORG  <  I  )  )  .GT.  XP.AX )  XI1AX  =  ABS  (XOP.G  ( I )  ) 

V ( I ) =VORG ( I ) 

X ( I ) =XORG ( I ) 

24  XXSU!l=XXSUMi-XORG  (I  )  *XORG  (I ) 

INOISE=0 

IFIVAR.GE.  l.E-10)  IUOISE=l 
IFQNOISE.  EQ.1)SNR  =  XXS(JM/(VAR*NPT) 

IF  (I  NOISE.  EO.  1 )  SMRDB  =  1 0 . 0  *ALOGl  0  ( SNR ) 

I  F  ( INOISE.  EQ. 1) WRITE  <6,1044  )  SNR,  SNRDB 
I F ( I NOISE. EO. 0 ) WRITE (6,1045) 

XV=0.5*XMAX/VMAX 
DO  22  K=1,MP1 
C  VORG(K)=XV*VORG(K) 

22  CONTINUE 

I F  (VAR .  G E .  1 . 0 E-  9 )  CALL  CORU  PT  ( XORG ,  X ,  V  AR ,  DPI ,  NAXPL ) 

C 

C  ADD  BIAS  TO  DATA  IF  IBIAS.NE.O 

C  (SET  B IASI  AND  B IAS 2  CONSTANTS  TO  ARTIFICIALLY  ADD  BIAS) 

I F ( IB  IAS, EQ. 0 ) GO  TO  6611 
BIAS1=0.0 
B IAS2=0 . 0 
DO  10  1  =  1,  MP1 
V  ( I )  =V  ( I )  +B  IAS  1 
10  X(I)=X(I)+BIAS2 

6611  CONTINUE 

C  IFdPR.GE.l.AND.  IPLT.  EO.  0 )  WRITE  (6 , 6993 ) 

C  IFdPR.GE.l.AND.  IPLT.  EQ.  0 )  1/RITE  (6,6994)  (X  (K )  , K=1 , IlPl ) 

C 

I F  (I PLT.  EQ .  0 ) GO  TO  6633 

C  SCALE  INPUT  WAVEFORM  FOR  REASONS  OF  PLOTTING  CONVENIENCE 

VSCAL=1 . 0 
DO  807  K=1,I1P1 
807  VORG(K)=VORG(K)*VSCAL 

IF  (ABS  (VSCAL-1 . 0 )  .GE.  1 .  E-6 )  1/RITE  (6 , 1047  )VSCAL 
IFdPR.GE.  1)  WRITE  (6, 1003) 

IF (IPLT. EQ . 1 . OR . IPLT. EQ. 3) 

+  CALL  PLOTIT  (DATA  ,  2 , [1P1 , 1  , MP1 ,  ISKIP,  MAXPL,  1,1.0) 

I F (IPLT. LE. 1 . OR. IPL1 . EQ. 0) GO  TO  6633 
C  LABEL  (18)  =CCHK 

6633  CONTINUE 

C  . 

C 

C  PERFORM  IDENTIFICATION  FROM  INPUT-OUTPUT  PAIR 

C  THEN  PREPARE  FILE  FOR  PLOTTING.  FILE  WILL  CONTAIN 
C  (NOISY  RESPONSE,  INPUT,  ORIGINAL  RESPONSE,  MODEL  RESPONSE) 

C 

CALL  GRAMI (X, V, MP1 , N, DELTA, Q, IZTS , GAMMA, XLAMDA, G, Z  ,  MAX, IREM, IDLY) 
CALL  ERROR ( XREC , VORG , GAMMA, DPI , N , XLAMDA , XORG ) 

IFdPLT.EQ.O)  GO  TO  6544 
IF (IPR.GE.l) WRITE (6,66) 

I F (IPR. GE. 1 ) WRITE (6,1004) 

IF (IPLT. EQ. l.OR. IPLT. EQ. 3) 

+CALL  PLOTIT (DATA 2 , 2 , MP1 , 1 , MP1 , ISKI P, MAXPL, 1,1.0) 

I F ( I PLT. EQ. 1 .OR. IPL2 . EQ . 0 ) GO  TO  6544 
C  LABEL ( 18 )=CCHK 

6998  FORMAT (2X, 'HERE  I  AM  IPLT=',l4) 

WRITE (9,6997) ( (DATA (K. 1 ) • DATA ( K . 2 ) . 
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+DATA2 (K,l) ,DATA2(K,2> ) ,K=1,MP1> 

6544  CONTINUE 
C 

C  WRITE (6,2) 

C  WRITE (6,1022) 

100  GO  TO  4321 
1234  CONTINUE 

C  . 

c 

c 

1  FOR11AT  !6I5,F10.0,6F5.0) 

2  FORUAT ( ' 1 ' ) 

3  FORUAT (/) 

66  FORUAT (//, IX, ’TRUE  RESPONSE  V  erSUS  RECONSTRUCTED  RESPONSE’,//) 

1000  FORUAT  (1H1,  SOX,  'STARTING  SIMULATION  ’, /20X,  •  SYSTEM  ORDER  =  ',15, 

1/ , 20X, ' M  +  1  =  ', 15,///, 20X, 'SAMPLING  INTERVAL  =  ' , FI 0 . 6 , // , 20X, 

2 ' IREM  =  ' , IS,/,20X, ' IB  IAS  »  ' , IS,/, 20X, ’ IDLY=  ’,15,/) 

1001  FORMAT (1215) 

1003  FORMAT  (//,5X,  '  INPUT  (  +  )  AND  OUTPUT  (*)  OF  THE  PLANT',/) 

1004  FORMAT  (//,  5X,  '  OUTPUT  (*)  AND  RECONSTRUCTION  (  +  )  ',/) 

1021  FORMAT (80A1) 

1022  FORMAT (/////////////////////////////////////////////////////////) 

1023  FORMAT (/ , IX, •*************************************************>, 

1  i **********************************  i ) 

1033  FORMAT ( 20 A4) 

1044  FORMAT <2X, 'SNR  =',F12.5,'  SNRDB= ' , F10 . 2 , / ) 

1045  FORMAT (2X, 'SNR  “INFINITY  (MO  NOISE)') 

1047  FORMAT (28, 'INPUT  WAS  SCALED, THEREFORE  DIVIDE  NUM  BY',/, 

+SX,F10.4) 

6994  FORMAT (2X, 10F10 . 4 ) 

6993  FORMAT ( 2 X, 'RESPONSE  DATA') 

6995  FORMAT (7F10.0) 

6996  FORMAT (2X,10F8.0) 

6997  FORMAT (6 (Gil. 4, IX) ) 

C  CALL  PICSIZ (0.0, 0.0) 

998  CONTINUE 

END  FILE  9 
STOP 
C 

C  DEFINITION  OF  PARAMETERS  USED  IN  THE  SIMULATION  OF  A 

C  LINEAR  DYNAMIC  SYSTEM 

C 

C  X  IS  THE  30RRUPTED  OUTPUT  SEQUENCE 

C  V  IS  THE  CORRUPTED  INPUT  SEQUENCE 

C  GAMMA  IS  THE  COEFFICIENT  VECTOR 

C 

C  MAX  =  ACTUAL  DIMENSION  SIZE  OF  2-DIM  ARRAYS  IN  THE  DIMENSION 

C  STATEMENT 

C  N  =  ORDER  OF  SYSTEM 

C  THE  MAXIMUM  VALUE  OF  N  IS  MAX/2-1 

C  MP1  =  M+l,  THE  TOTAL  NUMBER  OF  SAMPLED  POINTS  IN  EACH  SEQUENCE 

C 

C  RHO  =  EXPECTATION (  W(K)*Q(K)  ) 

C 

C  DELTA  IS  THE  SAMPLING  INTERVAL 

C 

C  IGRM  =  1  GRAMI  IS  PERFORMED 

C 

C  IPLT=0  NO  PLOTS 

C  I PLT=1  PLOTS  ONLY  WITH  PRINTER 

C 

C  IBIAS  -0  NO  BIAS  IS  ASSUMED  PRESENT  ON  INPUT-OUTPUT  DATA 

C  IBIAS  =1  SMALL  VALUES  OF  INPUT-OUTPUT  BIAS  ARE  ASSUMED  PRESENT 

C  ON  THE  DATA. 

C 
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subroutine  builda (A,Q, DEL, N, I1AX) 

REAL *8  A(I!AX,1)  ,Q(1)  ,DEL(1)  ,  PROD 

NP1=N+1 

NPNP2=N+N+2 

A ( 1 , NP1 ) =1 . ODOO 

PROD=1.0DOO 

D031 2I<=1 ,  N 

I=NP1-K 

PROD= PROD/DEL (I ) 

A (1> I) ’PROD 

312  A  (f\+l ,  NPl )  =0.  ODOO 
D0313  1  =  2 ,  NPl 
D0313K=1,N 

J=NP1-K  ,  % 

313  A(I,J)=(A(I, J+1)-Q(J) *A(I-1#J+1>)/DEL(J) 
D0314 1=1 , NPl 

D0314J=1 , NPl 
A(I, J+NP1)=0.0D00 
A  (I+NP1, J) =0. ODOO 

314  A(I+NP1,  J+IIP1)=A(X,  J! 

C  WRITE (6,1005) 

1005  FORfiATdX,  ’A-MATRIX') 

C  CALL  PR NAT ( A, NPNP2 , NPNP2 , MAX) 

RETURN 

END 

SUBROUTINE  BUILDZ  <Z , R, QI , NPl , NPT, NDIM) 

c  - 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  Z (NDIM, 1) ,R(1) 

COMMON  /DA2/IPR 
0=01 
N=NP1-1 
NPNP2=NPl+NPl 
R (1) =1. 0 
DO  1  1=1 ,  NPNP2 
IF(I.GE.2)R(I)=R(I-1) 

DO  1  J  =  1 ,  NPNP2 

1  Z<I,J)=0.0 
DO  2  K=1 ,  NPT 
NPTK=NPT+1-K 
DO  3  J=1,MP1 
DO  3  1= J ,  NPl 

3  Z(I,J)=Z(I,J)+R(I)*R(J)*HPTK 

R  (1 ) =0.0 
DO  4  1  =  1,  N 

4  R(I+1'=Q*R(I+1)+R(I) 

2  CONTINUE 

DO  174  1=1 , NPNP2 
DO  168  J= I , NPNP2 
168  Z(I,J)=Z(J,I) 

IF ( IPR.GE, 2) WRITE (6 , 220 ) <Z(I,J) ,J=1,NPNP2) 
174  CONTINUE 

220  FORMAT (2X, 5 (2X,G17 . 6) ) 

RETURN 

END 

SUBROUTINE  CORU PT  (F ,  X,  VAR , HPT,  NDIM) 

c 

c  CALLS  GGNML  (VIA  IMSL) 

C  ADDS  NOISE 

C  - - - - 

DIMENSION  F(l) ,X(1) ,R(1Q24) 

COMMON  /DA2/IPR 

DOUBLE  PRECISION  DT , DSEED 

FBAR=0  • 

EBAR=0. 

FESUM=0, 
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EESUtl=0. 

IF (IPR.GE.l) WRITE (6, 489) VAR 
C 

C  GENERATE  STANDARD  NORMAL  RV ,  THEN  CORRECT  STD  DEV 
C 

SIGMA=SQRT(VAR) 

DSEED= 123457. DO 

CALL  GGHt'.L  (DSEED,  NPT,  R) 

DO  26  K=1,NPT 
R(K)=SIGI1A*R(K) 

26  X(K)=F(K)+R(K) 

C 

DO  211  K=1,NPT 
FBAR=FBAR  +  F  (K) 

EBAR=EBAR+R(K) 

EESUM=EESUM+R(K) *R (K) 

FESl’M=FESUM+2.0*F (K) *R (K) 

211  CONTINUE 

FBAR  =  FBAR/NPT 
EBAK=EBAR/NPT 

IFdPR.GE.  1)  WRITE  (6,482)  FBAR, EBAR, FESUM, EESUM 
IF ( I  PR . LE. 2 )  GO  TO  411 
WRITE (6,8) 

WRITE (6,110) (X (K) , K=1 , NPT) 

WRITE (6, 18) 

WRITE  (6,115)  (R  ( K)  ,  K=1 ,  NPT) 

WRITE (6,1) 

411  CONTINUE 

999  CONTINUE 

C  FORMAT  STATEMENTS 

C 

8  FORMAT ( 10X, '  CORRU PTED  RESPONSE ' ) 

18  FORMAT  (1  OX,  '  ADDITIVE  NOISE  '  ) 

210  FORMAT ( 2X, 5 ( 2X, G14 . 7 ) ) 

110  FORMAT (20 (IX, F5.2) ) 

115  FORMAT (IX, 20 (IX, F5.3) ) 

482  FORMAT  (2X,  5HFB AR=  ,  Ell .4 ,6H  EBAR= ,  El  1 . 4 , 511  FE2  = ,  Ell .  4 , 4H  EE=,E11.4) 
489  FORMAT  (2X,  'VARIANCE  OF  NOISE=  '  ,  El  2 . 4  ) 

1  FORMAT!/) 

RETURN 

END 

SUBROUTINE  ERROR ( XREC,V , GAMMA , MP1 , N, XLAMDA, XORG , FDBACK , IDLY) 
DIMENSION  XREC(l) ,V(1) , XORG < 1 ) 

DIMENSION  WV  (  20  ) 

REAL  *8  GAMMA  (1)  ,  XLAMDA  (1)  ,  AVGW,  SUMV2  ,  AVGQ 

INTEGER  FDBACK 

NPNP2=N+N+2 

CALL  RESPOM (XREC,V, N, GAMMA, XLAMDA, MP1 , FDBACK) 

AVGW=0. 0D00 
SUMV2  =  0.0D0 
D026 1=1, MP1 

SUMV2  =  SUMV2  +  X0RG  (I)  *XORG  (I) 

AVGQ=XORG ( I ) -XREC ( I > 

26  AVGW'=AVGW+AVGQ*AVGQ 
AV  GW*  AV  GW/  SU  MV  2 
AVGQ=DSQRT (AVGW) 

AVGW=100.0*AVGW 

AVGQ=100.0*AVGQ 

WRITE  (6,27)  SUMV2,  AVGW,  AVGQ 

C27  FORMAT (IX, ' PER  CENT  MEAD  POWER  ERROR  OF  RECONSTRUCTION ', F8 . 3 ,/// , 
C  11X, ' PER  CENT  OF  SQUARE  ROOT  OF  POWER  ERROR  OF  RECONSTRU CT. 1 , F8 . 3 ) 

27  FORMAT ( IX,  ' SS  RESPONSE= ' , Ell . 4 , '  %  POWER  ERROR= ' , El  1 . 4 , / , 

+'  %  RMS  ERROR = ' , El 1 . 4 ) 

RETURN 

END 

SUBROUTINE  FILLV  (V,  NPT,  INPUT,  NPUL) 
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FILLS  THE  ARRAY  FOR  INPUT  ACCORDING  TO  INPUT  OPTION  DESIGNATED 
INPUT  =  0:  ZERO  INPUT 

1:  IMPULSE  AT  K  =  1 
2:  STEP  FN. 

3:  SQPULSE;  4:  DOUBLET,  PUL  WIDTH  =  NPUL 
5:  TR+  (  );  6:  TR+- (  ),  PUL  WIDTH  =  NPUL 

7:  SQUARE  HAVE,  PERIODS  PUL 

3:  SO.W  BURST  (P=NCY/l.’PUL,i;iDTH  =  NPUL)  ,  EXP  (TC=NPT/10) 

9:  EXP  (T.C.=NPUL) 

10:  IMP.  TRAIN'  <P=UPUL) 

11:  COS  OSC(  )  12:  SIN  OSC  T.C.=NPT/5,  P=NPUL 

13:  RANDOM.  S1GM.A=  .  1  *!!PUL 

14:  COS  BURST;  15:  SIN  B.  F  =  NCY/NPUL,  WIDTH=NPUL 
16:  COS  CHIRP  17:  SIN  C.  I  IN  IT  =  NCY/NPUL,  WIDTH  =  NPUL 

DIMENSION  V  < 1 ) 

COMMON  /WAVE/NCY 
DOUBLE  PRECISION  DSEED 
PI2=6. 28318530 
DO  90  1  =  1, NPT 
90  V(I)=0.0 

IF  ( INPUT.  EQ.  0)  GO  TO  999 

GO  TO  (1,2,3,3,5,5,7,8,9,10,11,12,13,14,15,16,17) , INPUT 

1  V ( 1 ) =1 . 0 

DO  101  1  =  2, NPT 

101  V(I)=0.0 
GO  TO  999 

2  CONTINUE 

I F  (NPUL.  EQ.  0)  NPUL=10 
AMP=0 . 1  *M PUL 
DO  102  1  =  1, NPT 

102  V ( I ) =AMP 
GO  TO  999 

3  CONTINUE 
Ml=NPUL/2+l 

DO  103  1=1, NPUL 
V  < I ) =1 . 0 

103  I F  (INPUT.  EQ.  4  .  AND.  I . GE.  N1 ) V  ( I )  =-l .  0 
GO  TO  999 

5  CONTINUE 

IF  (INPUT.  EQ.6)NPUL  =  NPUL/2 

Nl=NPUL/2+l 

N2=NPUL+1 

IF  (INPUT.  EQ.  6)  N2=NPUL+NPUL/2+l 

N3=NPUL*2+1 

NPUL2=NPUL/2 

DO  150  1=1, N1 

150  V ( I ) = FLOAT (I-1)/NPUL2 
DO  151  I=N1,N2 

151  V(I)=  2 . 0- FLOAT  ( 1-1 )  /NPUL2 
IF(INPUT.EQ.5) GO  TO  999 

DO  152  I=N2 , N3 

152  V ( I ) =-4 . 0+FLOAT ( 1-1 ) /MPUL2 
GO  TO  999 

7  DO  107  1=1, HPT 
V(I)=1.0 
TPUL=I/NPUL 

I F ( I-TPUL*NPUL . GE. NPUL/2 )  V(I)=-1.0 
107  CONTINUE 

GO  TO  999 

8  NP=NPUL/NCY 
MPD2=NP/2 

DO  55  1  =  1, NPUL 
V(I)=1.0 
IL=HOD ( I , NP) 

IF ( IL. GT. NPD2 ) V (I>=— 1.0 
55  CONTINUE 
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n  n  n  n  n  n 


108 

9 

109 

10 
no 

11 

in 

12 
c 

c 

112 

13 

113 

14 

114 

15 

115 

16 

116 

17 

999 


SN=S IGN  (1.0, V  ( N  PU  L )  ) 

TOO  .  2*  (HPT-NPUL) 

DO108I  =  NPUL,  t!PT 

ARC 1 = FLOAT  (NPUL  - 1 )  /TC 

V (I >=SN*EXP(ARG1) 

GO  TO  999 
DO  109  1=1, NPT 
ARG2=-FLOAT  ( I  )  /FLOAT  (NPUL) 

V ( 1 ) =  EXP (ARG2 ) 

GO  TO  999 

DO  110  1  =  1 ,  [IPT,  NPUL 
V(I)=1.0 
GO  TO  999 
DO  111  1  =  1 ,  !1PT 

AKG4=- 5 . 0  ‘FLOAT  (I )  /FLOAT  (NPT) 

ARG5=6 . 2832*FLOAT  (I )  /FLOAT  (NPUL) 

V  ( I )  =  EXP  ( ARG4  )  *COS  (ARG5) 

CONTINUE 

GO  TO  999 

DO  112  1=1, HPT 

AEG  3  = -FLOAT  ( I )  /FLOAT  (NPUL) 

ARG4=-5 . 0* FLOAT ( X ) /FLOAT (NPT) 

ARG5  =  6 . 283  2 ‘FLOAT  (I )  /  FLOAT  (NPUL) 

V  (I )  =EXP  ( ARG4  )  *S  It)  (ARG5) 

V (I ) =EXP ( ARG3 ) +EXP(AEG4) *SIN (ARG5) 

CONTINUE 

GO  TO  999 

DSEED= 789457 .DO 

CALL  GGHML (DSEED, NPT, V) 

IF  (NPUL.  EQ.  0  )  NPUL=  10 
DO  113  1=1, NPT 

V  (I )  =0. 1*V  (I )  *NPUL 
GO  TO  999 
CONTINUE 
NP=NR)L/NCY 

DO  114  1=1, NPUL 

V(I)=COS (FLOAT (1-1) *PI2/HP) 

GO  TO  999 

CONTINUE 

NP=NPUL/KCY 

DO  115  1=1, NPUL 

V (I ) =SIN (FLOAT (1-1) *PI2/NP) 

GO  TO  999 
CONTINUE 
WO=PI2*NCY/NPUL 
DO  116  1  =  1, NPUL 
TI2  =  2.0*  (I-U/NPUL 
WL= (1. 0+TI2**2) *WO 

V  ( I ) =SIN (FLOAT (1-1 ) *WW) 

GO  TO  999 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE  FIX  (G,  P,  C,  D,  X,  N,  HC,  SIG,NDIM,  IFIX) 

IMPLICIT  REAL*8  (A-H,0-Z) 

CORR^Vo?sYEr.ATRIxfICY  SIG  ‘ASSUf,E  "HITE  N°ISE> 

P  DENOTES  NOISE  MATRIX  FOR  UNIT  NOISE 

NC  IS  THE  NONZERO  SUBMATRIX  OF  P  =COV  OF  NOISE 

DIMENSION  G (NDIM, 1 ) ,P(NDIM,1. ,C(NDIM,1> , D (NDIM, 1 ) , X ( 1) 
DOUBLE  PRECISION  SUMDET, DT 
COMMON  /DA2/IPR 
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COMMON  /DA3/DFAC 
SI=SIG 

IF(IFIX.EQ. 0)GO  TO  51 
GDP«G(1,1)/P(1,1) 

JCT=0 
SIG  =  0.0 

210  FORMAT (IX, 5G12.5) 

DO  211  1=1,  N 

211  IFfIPR.GE.2)  WRITE  (6,210)  (G(1,J)  ,J=l,tl> 

3  JCT=JCT+1 

SUMDET=0.0 

CALL  GKKDCT  (G,  D,  GDET,  X,  N,  tJDItl,  0} 

DO  212  1=1, N 

212  IF(IPR.GE.2)WRITE<6,210)  (D  ( I ,  J  )  ,  J  =  1 ,  N) 

I F ( JCT. EQ. 1 ) DETG=GDET 

DO  7  1=1, NC 
DO  7  J=1,NC 

7  SUMDET=SUI1DET+D(I,J>  *P(I,J) 

I  F  (SU  tlDET.  LT.O.O.AIiD.  IFIX.NE.2)  GO  TO  11 
S  1  =  1 .  0D0 /SUMDET 

WRITE (6,32) JCT, ICT.GDET, SUMDET, SI 
IF'SI/GDP.CT.O.DWRITE  (6,31) 

ICT=0 

51  CONTINUE 

DO  9  1=1, N 
DO  9  J  =  1,N 

C(I,J)=G(I,J)-SI*P(I,J) 

9  CONTI HUE 

I F ( I  FIX.  EQ. 0 ) GO  TO  11 

CALL  GKRDCT  (C,  D,  CDET,  X,  N,  NDII1,  0) 

I F (CDET. LT. 0 . 0 . OR . CDET. GT. GDET) ICT=ICT+1 
IF  (ICT.GT.  5)  GO  TO  10 
IF ( ICT.GT. 0)SI=SI/2.0 
IF (ICT.GT. 0 ) GO  TO  51 
IF(JCT.GE.  S)GO  TO  11 

10  THR=DETG*DFAC 

I F ( J  CT . EQ. 1) WRITE (6,33) DETG, DFAC, THR 
SIG=SIG+SI 

I F (CDET. LE . THR ) GO  TO  22 
DO  23  11  =  1, N 
DO  23  JJ=1,N 
23  G(II,JJ)=C(II,JJ) 

22  CONTINUE 

WRITE (6,34) JCT, ICT, GDET, SI , CDET 
I F (CDET. GT. THR ) GO  TO  3 

31  FORMAT  <2X,  '  NOISE  VAR  EXCESSIVE,  SIG/GDP. GT.  0 . 1 1 ) 

32  FORMAT (IX,  'J, I  GDET, SUMDET, SI 2 12 , 4 Ell . 2 ) 

33  FORMAT (IX, 'GDET, DFAC, THR: ' ,6E11.2> 

34  FORMAT (IX, 'J, I  S I , CDET ' , 2 12 , 4E11 . 2) 

11  CONTINUE 

DO  25  1=1, N 

IF ( I  PR. GE. 2) WRITE (6, 210)  (C ( I , J ) , J=1 , N) 

DO  25  J=1,N 
25  G(I,J)=C(I,J) 

RETURN 

END 

SUBROUTINE  GKRDCT (X, Y, DET, XLAHDA, N, MAX, IOPT' 
REAL*8  X  (MAX,  1 )  ,Y(MAX,1)  ,  A,  B ,  C,  D,  E,  DET,  XLAMDa  i.  ) 
REAL  *8  SCAL  (20)  ,RSC(20)  ,DT,SC,  AD,  BD 
INTEGER  NUI1  (2,20) 

COMMON  /DA2/IPR 
COMMON  /GKRD/IGKR 
COMMON  /GKRD2/DT 

C  IGKR=0  USE  IS  MADE  OF  THE  FIRST  ROW  OF  ADJOINT 

C  1  DIAGONAL (NEGATIVE  ENTRIES  SET  TO  ZERO) 

C  2  ABSOLUTE  VALUE  OF  DIAGONAL 
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c 

C  IOPT  =  0:  RETURN  X  INVERSE  IN  Y 

c  i:  Calculat#ARAMETER  VECT0R 

I F (N . NE. 1 )GO  TO  3 
YU,1)=1.0/XU,1> 

DET=Xd,  1) 

GO  TO  61 
3  CONTINUE 

C 

C  SCALE 

C 

ISCL=1 

IFdSCL.EQ.O) WRITE (6,804) 

804  F0R!1AT(2X,  '** ‘WARNING ! !  SCALING  IN  GKRDCT  DISABLED**') 

DO  11  1  =  1, N 
SCAL(I)=1.0D0 

IFdSCL.GE.l.AND.Xd,  I)  .GT.  0 . 1 E-20  )  SCAL  ( I )  =DSQRT  (X  ( I ,  I)  ) 
11  RSC (I ) =1 . 0/SCAL (I ) 

DO  6  1  =  1,  N 
DO  6  J=1,N 

6  Y(J,I>=X(J,I)*PSC(I)*RSC(J) 

C  CALL  PRMAT  (Y,  N, N, MAX) 

A-1.0D0 
DO  43  1=1, N 
B=0.0D0 
L=  I 
11=1 
C 

C  FIND  LARGEST  ENTRY  A  (L,  (1)  IN  LOWER  DIAGONAL  SUBNATRIX 

C 

DO  18  J= I , N 
DO  18  K= I , N 

IF  (DABS  <Y  (tl,  J 1  )  .  LE.B)  GO  TO  18 
B=DABS ( Y (K,J) ) 

L=K 

H=J 

18  CONTINUE 
C 

C  INTERCHANGE  ROWS 

C 

IF1L.EQ.  I ) GO  TO  24 
DO  23  J=1,N 
C=Y (L, J) 

Y(L,J)=Y(I,J) 

23  Y ( I , J) =C 
C 

C  INTERCHANGE  COLUMNS 

C 

24  IFU1.EQ.  DGO  TO  29 
DO  28  J= 1 , N 

C= Y ( J , M) 

Y  (J,M)=Y(J,  I) 

28  Y ( J , I ) =C 
C 

C  BEGIN  SWEEP  COLUMNS  TO  THE  RIGHT 

C  ARRAYS  NUM ( 1 , - )  ,NUM(2,.)  KEEP  RECORD 

C  OF  ROW  AND  COLUMN  INTERCHANGES 

C 

29  NUM  < 1 , 1 ) =L 
NUM (2 , 1 )  =M 
B=Y  (1,1) 

Y(I,  I)=A 

DO  42  J=1 , N 

IF  (J.  EQ.  DGO  TO  42 

C=-Y (I , J ) 
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Y(I,J)=O.ODO 
DO  41  K  =  1,N 
D=Y (K, I) »C 
E=Y(K, J) *D+D 

C  I F ( DAB  S ( E ) .  LT. 1 . 0D-10 *DABS <D) )E=0.0D0 

41  Y(K,J)=E/A 

42  CONTINUE 

43  A=B 
C 

C  RESTORE  COLUMNS 

C 

DO  58  1=2, N 
J-N+l-I 
K=NUM(2,  J) 

IF  (K. EQ. J ) GO  TO  52 
DO  51  L=1 ,  N 
C=Y(K,L) 

Y(K,L)=Y(J,L) 

51  Y  ( J ,  L)  =C 

52  K=NUM ( 1 , J ) 

C 

C  RESTORE  ROWS 

C 

IF(IC.EQ.J)GO  TO  58 
DO  57  L=1 , N 
C=Y (L, K) 

Y(L,IO=Y(L,J) 

57  Y (L, J ) =C 

58  CONTINUE 
DET=A 

DO  59  1=1, N 

59  DET=DET*SCAL  (I )  *SCAL  ( I ) 

I F  < IPR. EQ. 3 ) WRITE (6,337) DET, A, (RSC(I) , 1=1, N) 

IF  (IOPT.  EQ.  1)  GO  TO  61 
DO  60  1=1, N 
DO  60  J=1,N 

60  Y(I,J)=Y(I,J)*RSC(I)*RSC(J)/A 

61  CONTINUE 

C  WRITE (6, 703) 

703  FORMAT (2X, 'GKRDCTj  PROCESSED  MATRIX  Y’> 

C  IF (IPR. EQ. 1) CALL  PRMAT ( Y, N, N, MAX) 

IF  ( IOPT.  NE.  1 )  GO  TO  1000 
IF(Y(1,1) .LT.0.0D0)  GO  TO  1000 
C 

DO  806  1=1, N 

806  IF  (IPR.GE. 3) WRITE (6,803) (Y(I,J) ,J=1,N> 

SC=1 . ODO 
DO  200  1=2, N 
SC=SC*DT 

RSC ( I ) =RSC ( I ) /RSC ( 1 ) 

A=Y(I,I) 

803  FORMAT (2X, 6G12. 5) 

XLAMDA ( I >  =RSC (I ) *DSQRT ( A/Y (1,1)) *DSIGN ( 1 . ODO , Y ( 1 , 1 ) ) 

200  CONTINUE 

XLAMDA ( 1 ) = 1 . ODO 

IF (IPR. EQ. 1 ) WRITE (6,106) (XLAMDA (IP) , IP=1 ,N> 

106  FORMAT (5X, ' SYNTHETIC  PARAMETER  V ECTOR ’, 10G12 . 5 ) 

1000  CONTINUE 

337  FORMAT (IX, 'DET, A, RSC (I) :  ',7E11.2) 

339  FORMAT (IX, 'ERROR.  RATIO  OF  11=', II,'  JJ=',I1,'  COFACNEG' 
RETURN 
END 

SUBROUTINE  GRAMI (X,V, MP1 , N, DELTA, QSAV , KOPT, GAMMA, XLAMDA, G,Z 
1IREH, IDLY) 

C  THIS  SUBROUTINE  PERFORMS  THE  GRAMI  TECHNIQUE 
C 


i 


) 

,  MAX, 
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DIMENSION  XU) ,V(I) ,G (MAX, 1 ) , Z (MAX, 1 ) .GAMMA (1 ) , XLAMDA ( 1 ) ,Q(20) , 
1DEL ( 07 ) 

DIMENSION  GAIH25) 

DOUBLE  PRECISION  GAM 

DOUBLE  PRECISION  G, Z, GAMMA, XLAMDA, DELTA, DEL, PROD,Q,QSAV 
DIMENSION  GDUI1  (20,20)  ,GEST(20,20)  ,GGG(20,20> 

REAL'S  GDUI!,GEST,GGG 
REAL *8  VAP.Q.VARW,  FAC 
REAL'S  SCALE (20) , SCAL 
COMMON  /GKRD/IGKR 
COMMON  /GKRD2/DT 
COMMON  /DA1/NN, IFIX 
COMMON  /DA2/IPR 
COMMON  /DA3/DFAC 
COMMON  /B IAS/ IB  IAS 
C 

DT=DELTA 

I F ( I  PR. GE. 1) WRITE (6,1000) 

1000  FORMAT (1H1.20X, ’THE  GRAM  I  TECHNIQUE’) 

C  JOPT  =  0  IF  DIRECT  TRANSMISSION  IS  ASSUMMED 

JOPT=0 

IF(IREM.NE.O) JOPT=l 
C 

C  DEL  IS  THE  NUMERATOR  OF  THE  KNOWN  FIRST  ORDER  DIGITAL  FILTERS 

D019I=1,N 
DEL ( I ) =1 , 0D00 
19  Q(I)=QSAV 

IF(IPR.GE.O) WRITE (6,2020) QSAV 
2020  FORMAT (3 OX, ’Q  PARAMETERS:  Q=’,F8.3) 

C  CALL  PRVEC(Q,N) 

NP1-N+1 
NP2=N+2 
NPUP1-N+N+1 
NPNP2=N+N+2 
NR=NPl-IREM 
HP1PIR=NP1+IREM 
DO  12  1=1, MAX 
DO  12  J=1 , MAX 
12  Z (I,J) =0.0D00 

VARO=0.O 
VARW=0.0 
DO  300  1=1, MP1 
VARW=VARV7+V  (I)  *V  (I) 

300  VARQ=VARQ+X(I)*X(I) 

VARQ=DSQRT  (VARQ/MP1 ) 

VARW=DSQRT  (VARW/MP1 ) 

C 

C 

C  CALCULATING  THE  G  MATRIX 

IF ( IB IAS. EQ. 0 ) GO  TO  11 
HPNP2=N+N+2+l 
NR-NPl-IREM+l 
11  CONTINUE 

C 

DO10I=l ,NPNP2 
GAM (I) =0.0 
GAMMA ( I ) =0 . 0D00 
DO10J=  I ,  NPNP2 
10  G ( I , J) =0 . 0D00 

GAM (1) =1.0 
DO50K=l , MP1 
IF (K- IDLY  >25,25,24 
25  GAMMA (NP2)=0.0D00 

GO  TO  26 
24  FAC= 1.0 

GAMMA (NP2) =V (K-IDLY) /FAC 
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FAC-1.0 

GAMMA (1)=X(K> /FAC 
26  CONTINUE 
DO30 1=1 , N 

GAM ( 1+1 ) =GAM (1+1 ) *Q ( I ) +GAM (1 ) ‘DEL ( I ) 

GAMMA ( 1+1 ) =GAMMA ( I ) *DEL ( I ) +GAMMA ( 1+1 ) *0 (1 > 

30  GAMMA ( I+HP2 ) =GAHMA ( 1 +HP1 ) ‘DEL ( I ) +GAMMA ( I+HP2 )  *Q  (I  ) 

IF(IBXAS.E0.  1)  GAMMA  (NPHP2  >  =GAM  ( 1  +  1  > 

DO  40  1  =  1  /  NPNP2 
DO  40  J= I # NPNP2 

40  G(I,J)=G(I,J)  +GAI1MA  { I )  ‘GAMMA  ( J ) 

C  WRITE (6, 1025) K, (GAMMA (I) ,  I=1,NPNP2) 

1025  FORMAT (2X,10F8.3) 

50  CONTINUE 

C 

C  PERFORM  SCALING  ON  G-MATRIX 

DO  735  1=1 , NPNP2 
SCALE ( I ) =DSQRT (G ( I , I ) ) 

SCALE  (I )  =1 . 0 

735  CONTINUE 

C  WRITE (6 , 1008) 

1008  FORMAT (10X, 'SCALE  VECTOR') 

C  CALL  PRVEC (SCALE , NPNP2 ) 

DO  736  1=1 , NPNP2 
DO  736  J=1 , NPNP2 

736  G(I,J)=G(I,J)/ (SCALE ( I ) ‘SCALE ( J ) ) 

1023  CONTINUE 

C  WRITE (6, 1009) 

1009  FORMAT  (10X,  ' - G  MATRIX - *) 

C  DO  56  1=1 , NPNP2 

C56  IF (IPR.GE. 2) WRITE (6,3) (G(J,I) ,J=1,I> 

3  FORMAT (IX,I0D13.5) 

42  CONTINUE 

DO60I=2,NPNP2 

K=I-1 

DO60 J=1 , K 

60  G(I,J)=G(J,I) 

C 

c 

C  NOISE  ESTIMATION 

IF  (IFIX.  LE.  0) GO  TO  878 

CALL  BU ILDZ ( Z , XLAMDA , QSAV , NP1 , MP1 , MAX ) 

CALL  FIX (G, Z ,GEST, GDUM, XLAMDA, NPHP2 , NP1 , SIG2 ,MAX, 1 ) 
878  CONTINUE 

C 

c 

IF (JOPT) 70,90,70 
70  DO80J=l , NPNP2 

DO80I=l , NR 

G (NP1  +  I, J) =G (NP1PIR+I,  J) 

80  CONTINUE 

NPNP2=NPNP2-IREM 
D085J=1 ,NPNP2 

SCALE (NP1+J ) “SCALE (NP1PIR+J) 

D085 1=1 , NR 

G ( J , NP1+I ) =G ( J » NP1PIR+I ) 

85  CONTINUE 

90  CONTINUE 

DO  55  1=1 , NPNP2 

55  IF (IPR.GE. 2 (WRITE (6, 3) (G ( J , I ) , 0= 1 , I ) 

CALL  GKRDCT  (G,  Z ,  DET,  XLAMDA ,  NPNP2 ,  MAX,  1 ) 

C  DE-SCALE  SYNTHETIC  COEFFICIENT  VECTOR,  XLAMDA 
DO  741  I=1,NPNP2 

741  XLAMDA (I) -XLAMDA (I) /SCALE (I) 

I F ( IB  IAS . EQ. 0 ) GO  TO  43 
NPNP2-NPNP2-1 


A0 


l 


NR=NR-1 

43  CONTINUE 

XMEAN=XLAMDA (NPNP2+1) 

IF(JOPT) 120,130,120 
120  NPt!P2=NPNP2  +  IREM 

I F  ( I  PR.  GE.  2 )  WRITE  (6,210)  (XLAMDA  (I)  ,  I=1,NPNP2) 

210  FORMAT (IX, 5G12 . 5 ) 

D0122 1=1, NR 

122  XL AM DA (NPNP2-I  +  1 ) =XLAMDA (NP2 VNR- I > 

DOl 23  1=1,  IREtl 

123  XLAMDA (NP1+I ) =0 . 0D00 
130  CONTINUE 

FAC=1 . 0 

DO301I=NP2,NPNP2 
301  XLAMDA  ( I )  =XLAMDA  ( I  >  *FAC 

IFdPR.GE.  3) WRITE  (6,1001) 

1001  FORMAT (10X, 'THE  SYNTHETIC  COEFFICIENT  VECTOR,  XLAMDA,  IS') 
IFdPR.GE.  3)  CALL  PRV  EC  (XLAMDA,  NPNP2  ) 

DO  150  1=1, NPNP2 
150  GAMMA ( I ) =XL AMDA ( I ) 

C 

C  GENERATING  GAMMA  FROM  XLAMDA 

CALL  BUILDA (G, Q, DEL, N, MAX) 

DO160I=l,NPNP2 
GAMMA (I)=O.ODOO 
DO160J=l , MPNP2 

160  GAMMA ( I ) =GAMMA (I )  +G  (I ,  J ) ‘XLAMDA ( J ) 

165  CONTINUE 

D0200 1  =  2 ,  NPNP2 

200  GAMMA  ( I )  =GAM11A  (I)  /GAMMA!  1) 

GAMMA ( 1 ) =1 . 0D00 

IF  (IDLY.  EQ.  0)  GO  TO  172 

IDLY1=IDLY+1 

DO  170  I 1= IDLY1 , NP1 

I=NPNP2+1-II 

170  GAMMA ( I+IDLY) =GAMMA ( I ) 

DO  172  1=1, IDLY 
GAMMA (I+NP1 ) =0 . 0D00 
172  CONTINUE 
C 

C  CALCULATING  THE  EQUIVALENT  CONTINUOUS  DESCRIPTION 

C  CALL  PRVEC (GAMMA, NPNP2) 

CALL  IZTOS (GAMMA, N, DELTA, KOPT) 

IF (IPR.GE.l) WRITE (6,1003) 

1003  FORMAT (/// ,1X,100(1H-) , / , IX, 100 ( 1H-) ) 

RETURN 

END 

SUBROUTINE  IZTOS  (GAMMA,  N,  DELTA,  IZTS) 

C 

C  IZTOS  SEPARATES  THE  NUMERATOR  FROM  THE  DENOMINATOR  PARAMETERS 

C  IN  GAMMA 

C 

DIMENSION  GAMMA ( 1) ,X1 (20) ,X2 (20) 

DOUBLE  PRECISION  GAMMA, XI , X2 , DELTA 
COMMON  /DA1/NN, IFIX 
COMMON  /DA2/IPR 
NP1=N+1 

200  D03 1=1 ,NP1 

XI  (I)=GAMMA(I) 

3  X  2 ( I ) =-GAMMA (NP1  +  I ) 

CALL  ZTOS(Xl,X2,N, DELTA, IZTS) 

IZTS=IZTS+1 

IF(IZTS.EQ.4)  GO  TO  200 

RETURN 

END 

SUBROUTINE  POLCON  (C,  R2,  K,N) 
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non  nnnnnnnnnnn  nn  nnn  n  nn  nnn  nnn  n  n  n  n  nn  t*  tj  ►-  nnn 


a  POLYNOMIAL  construction  program  needed  for  ztos 

DIMENSION  C ( 1 ) ,R2(1) 

COMPLEX* 16  C,R2,C0MP 
REAL *8  DC (2) 

EQUIVALENCE  (COMP,  DC) 

NP1=N+1 
D010 1=2 ,  NP1 
)  R2 ( I ) =0 . 0D00 

P,2(1)=1.0D00 
D04 1=1 , U 
COMP=C ( I ) 

IF  <  I  .  EQ.K.OR.  (DC  (1)  .EQ.0.0D0.  A1)D.DC(2)  .  EQ.0.0D0)  )  GO  TO  4 

D02JJ=1 , I 

J=I-JJ+1 

R2  (J+l )  =R2  (J  +  l)  *C(I)+R2  (J) 

R2  (1 ) =R2 ( 1) *C (I ) 

CONTINUE 

RETURN 

END 

SUBROUTINE  POLRT  (XCOF,  COF,  [1,  ROOTR ,  ROOTI ,  IER  ) 

COMPUTES  THE  REAL  AND  COMPLEX  ROOTS  OF  A  REAL  POLYNOMIAL 
DESCRIPTION  OF  PARAMETERS 

XCOF  -VECTOR  OF  M+l  COEFFICIENTS  OF  THE  POLYNOMIAL 
ORDERED  FROM  SMALLEST  TO  LARGEST  POWER 
COF  -WORKING  VECTOR  OF  LENGTH  M+l 
M  -ORDER  OF  POLYNOMIAL 

ROOTR-RESULTANT  VECTOR  OF  LENGTH  M  CONTAINING  REAL  ROOTS 
OF  THE  POLYNOMIAL 

ROOTI- RESULTANT  VECTOR  OF  LENGTH  M  CONTAINING  THE 

CORRESPONDING  IMAGINARY  ROOTS  OF  THE  POLYNOMIAL 
IER  -ERROR  CODE  WHERE 
IER=0  NO  ERROR 
IER=1  M  LESS  THAN  ONE 
IER=2  M  GREATER  THAN  36 

IER=3  UNABLE  TO  DETERMINE  ROOT  WITH  500  INTERATIONS 
ON  5  STARTING  VALUES 
IER=4  HIGH  ORDER  COEFFICIENT  IS  ZERO 

DIMENSION  XCOF ( 1 ) ,COF(l) ,ROOTR(l> , ROOTI (1) 

DOUBLE  PRECISION  XO, YO, X, Y, XPR, YPR,UX,UY, V, YT, XT, U, XT2 , YT2 , SUMSQ, 

1  DX, DY, TEMP, ALPHA, XCOF, COF, ROOTR, ROOTI 

LIMITED  TO  36TH  ORDER  POLYNOMIAL  OR  LESS. 

FLOATING  POINT  OVERFLOW  MAY  OCCUR  FOR  HIGH  ORDER 
POLYNOMIALS  BUT  WILL  NOT  AFFECT  THE  ACCURACY  OF  THE  RESULTS. 

METHOD 

NEWTON-RAPHSON  ITERATIVE  TECHNIQUE.  THE  FINAL  ITERATIONS 
ON  EACH  ROOT  ARE  PERFORMED  USING  THE  ORIGINAL  POLYNOMIAL 
RATHER  THAN  THE  REDUCED  POLYNOMIAL  TO  AVOID  ACCUMULATED 
ERRORS  IN  THE  REDUCED  POLYNOMIAL. 


I FIT=0 

N=M 

IER=0 

IF (XCOF (N+l ) ) 10,25,10 
10  IF (N)  15,15,32 

SET  ERROR  CODE  TO  1 

15  IER-1 

20  TF(IER)2O0.201,200 
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200 

203 

201 


WRITE (6,203) IER 

FORMAT (IX, 'ERROR  CALLED  FROM  POLRT,  IER  = 
RETURN 

SET  ERROR  CODE  TO  4 


25  I£R=4 
GO  TO  20 

SET  ERROR  CODE  TO  2 

30  IER=2 
GO  TO  20 

32  IF  <  13— 36)  35,35,30 
35  NX=N 
NXX=N+1 
N2=l 

KJ1  =  H+l 
DO  40  L=1 ,  KJ1 
MT-KJ1-L+1 
40  COF (MT ) =XCOF (L) 

SET  INITIAL  VALUES 

45  XO=. 00500101 
YO=0. 01000101 

ZERO  INITIAL  VALUE  COUNTER 

IN=0 
50  X=XO 

INCREMENT  INITIAL  VALUES  AMD  COUNTER 

XO=~10.0*YO 

YO=-10.0*X 

SET  X  AND  Y  TO  CURRENT  VALUE 

X=XO 
Y=YO 
IN=IN+1 
GO  TO  59 
55  IFIT=1 
XPR=X 
YPR=Y 

EVALUATE  POLYNOMIAL  AND  DERIVATIVES 

59  ICT=0 

60  UX=0.0 
UY=0 . 0 

v  =0.0 

YT=0 . 0 
XT=1 . 0 
U=COF (N+l ) 

I F (U )  65,130,65 
65  DO  70  1  =  1, N 
L  =N-I+1 
TEMP=COF(L) 

XT2=X*XT-Y*YT 
YT2=X*YT+Y*XT 
U=U+TEMP*XT2 
V=V+TEMP*YT2 
FI  =  I 

UX=UX+FI*XT*TEMP 


M3) 


A3 


nnnn  t-  nnn  non- j 
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U Y«UY-FI *YT*7EHP 
XT“XT2 
70  YT= YT2 

SU  MSQ-U  X*U  X+U  Y  *U  Y 
IF(SUMSQ)  75,110,75 
75  DX=  <V*UY-U*UX) /SUMSQ 
X=X+DX 

DY  =  -  (U*UY4V»UX)  /SUMSQ 
Y= Y+DY 

8  IF  (DABS (DY)4DABS(DX)-1.0D-10>  100,80,80 

STEP  ITERATION  COUNTER 

80  ICT= ICT+1 

IF ( ICT-500 )  60,85,85 
85  IF(IFIT)100,90,100 
90  IF ( IN-5 >  50,95,95 

SET  ERROR  CODE  TO  3 

95  IER=3 
GO  TO  20 

100  DO  105  L=1,NXX 
(1T=KJ  1-L+l 
TEMP=XCOF (MT) 

XCOF(rtT)=COF(L) 

105  COF(L)=TEIlP 
ITEIIP=N 
N-NX 

NX=ITEMP 

I F ( I  FIT)  120,55,120 
110  IF(IFIT)  115,50,115 
115  X=XPR 
Y=YPR 
120  IFIT=0 

22  IF (DABS ( Y) -1 . 0D-8*DABS (X) >135,125,125 
125  ALPHA=X+X 

SUMSQ=X*X4Y*Y 
N=N-2 
GO  TO  140 
130  X=0.0 
NX=NX-1 
NXX=NXX-1 
135  Y=0 . 0 

SUMSQ=0.0 

ALPHA=X 

N=N-1 

140  COF  (2)  =COF  ( 2 )  +ALPHA*COF  (1 ) 

145  DO  150  L=2,N 

150  COF  (L+l )  =COF  (L+l )  4ALPHA*C0F  (L)  -SUMSQ*COF  (L-l ) 

155  ROOTI (N2 >  =Y 
ROOTR (N2 ) =X 
N2=N24l 

IF(SUMSQ)  160,165,160 
160  Y=-Y 

SUMSQ«0.0 
GO  TO  155 
165  IF(tl)  20,20,45 
END 

SUBROUTINE  PRCVEC(A,N) 

PRCVEC  PRINTS  A  DOUBLE  PRECISION  COMPLEX  VECTOR 

C  A  COMPLEX  NUMBER  OF  THE  FORM  A+JB  IS  PRINTED  (  A,  B  J) 

DIMENSION  A ( 1) 

COMPLEX*16  A 
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I 


IF(N.EQ.U)  GO  TO  100 
WRITE (6,920) 

WRITE (6,910) (A (I ) , 1=1 , N) 

910  FORMAT (1X,1H(,D22.15,1H,,D22.15,3H  J ) ) 

C  WRITE (6, 920) 

100  CONTINUE 
920  FORMAT ( 2 ( / ) ) 

RETURN 

END 

SUBROUTINE  PRHAT  (A,  N,  M,  UMAX) 

DOUBLE  PRECISION  A 
C 

C  THIS  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  DOUBLE  DIMENSIONED  ARRAY 
DIMENSION  A (NMAX, 1 ) 

WRITE (6,1) 

D02 1=1 , N 

2  WRITE (6,3) <A(I,J) ,J=1,M) 

3  FORMAT (1X,10D13.5) 

WRITE (6,1) 

WRITE (6,1) 

1  FORMAT (/) 

RETURN 

END 

SUBROUTINE  PRVEC(A,N) 

C 

C  THIS  SUBROUTINE  OUTPUTS  DOUBLE  PRECISION  SINGLE  DIMENSIONED  ARRAY 

DIMENSION  A(l) 

DOUBLE  PRECISION  A 
WRITE  (6,1)  (A(I), 1=1, N) 

1  FORMAT (1X,10D13*5) 

31  FORMAT (/) 

RETURN 

END 

SU B  ROU  TINE  RES  PON ( X , V , N , GAMMA , XLAMDA, MP1 > 

DIMENSION  X<  1)  ,V(1)  , GAMMA (1 )  , XLAMDA  (1) 

REAL *8  XSAV, GAMMA, XLAMDA 

NH1-H-1 

NP1=N+1 

NPNP1=N+N+1 

NPNP2=N+N+2 

DO  19  1=1 , NPNP1 

19  XLAMDA (I)=0.0D00 
XSAV=0 . 0D00 

C 

C 

DO  20  K=1 , MP1 
IF (N.  EQ.  DGO  TO  25 
DO  21  1=1, NM1 
J=NP1-I 

21  XLAMDA (J)=XLAMDA(J-1) 

25  CONTINUE 

DO  22  1=1, N 
J=NPNP2-I 

22  XLAMDA ( J ) =XLAMDA ( J-l ) 

XLAMDA (1) =XSAV 
XLAMDA (N PI )=V(K) 

XSAV=0 • 0D00 

DO  23  1=1  ,NPNP1 

23  XSAV= XSAV -GAMMA ( 1+1 ) ‘XLAMDA (I ) 

IF  (DABS  (XSAV)  .GE.  1 . 0D10  )  XSAV  =  0 . 0D00 

20  X (K ) =XSAV 

77  FORMAT (2X,10F7.2) 

78  FORMAT  (IX,/) 

RETURN 

END 

SUBROUTINE  ZTOS (B, A, N, DELTA, IZTS) 


COMMON  /DA2/IPR 
COMMON  /DELAY/ IDLY 


C 

C 

C  CONVERSION  OF  A  DISCRETE  TIME  SYSTEM  H(Z>  TO  A  CONTINUOUS  TIME  SYS 

C 

C 

C  H (  Z )  =  ( A  ( 1 )  +  A  (  2 )  *  Z  ETA  +....)/(l  +  B(2)*ZETA  +....) 

C  ZETA  =  1/Z 

C 

C  H(S)  =  (A<1>  +  A ( 2 ) *S  +  .  +  A (N+l ) *S**N ) /DENOM 

C 

C  DENOM=B ( 1 )  +  B ( 2) *S  +  .  +  B(N+1)*S**N 

C 

C  B( 1)  =  1  ALWAYS 

C 

C 

DIMENSION  B ( 20 ) ,A(20) ,TEHP(20) ,RR(20) ,RI (20) ,CR(20) ,CA(20) , 

+CAA ( 20 ) ,CA1 (20) ,CB(20) ,CF(20) ,CF1 (20) ,CG(20) 

COMPLEX*16  CA, CAA,  CA1 ,  CB ,  CR, CONI ,CON2,CONT, FAC, A1 , A2 , B1 ,B2 , AA1 , BB1 
ICG , CF1 , CF 

REAL *8  B, A, TEMP, RR,RI, DELTA 

CONT=0 . 0D00 

IORP=  IZTS 

NP1=N+1 

NNP1=NN+1 

IF(IPR.GE.l) WRITE (6,989)  HN 
989  FORMAT (I0X, 'NN  =  ',15! 

999  FORMAT (/) 

IF (IPR.GE.l) WRITE (6,999) 

WRITE (6,1000) 

1000  FORMAT ('  Z-DOMAIN  DENOMINATOR') 

CALL  PRVEC (B, NP1) 

WRITE (6,1001) 

1001  FORMAT ('  Z-DOMAIN  NUMERATOR') 

CALL  PRVEC (A, MPI ) 

IF(IZTS.EQ.O)  GO  TO  919 

IF ( IZTS. EQ. 1 )  GO  TO  200 
IF(IZTS.E0.2)  GO  TO  250 
IF (IZTS. EQ. 3)  GO  TO  200 
IF (IZTS. EQ. 4 )  GO  TO  250 
200  CONTINUE 
C 

c 

c  LOGARITHMIC  TRANSFORMATION 

C 

C 

C  WORK  ON  NUMERATOR 

C 

C 

NMD=NP1-  IDLY 
DO  14  1=1, NMD 

14  A(I) =A (1+ IDLY) 

A1=0 . 0D0 

B1=0 . 0D0 

DO  301=1, NP1 

IF (I . LE. NNP1 ) A1=A1+A ( I ) 

30  B1=B1+B ( I ) 

I F (NN • EQ. 0 ) GO  TO  469 

CALL  POLRT (A, TEMP, NN, RR, RI, IER ) 

D015  1=1, HN 

15  CA ( I ) =DCMPLX (RR ( I ) , RI ( I ) ) 

DO  7  1=1, NN 

7  CA ( I ) = ( +1 . 0 /DELTA) *CDLOG (CA ( I ) ) 

IF (NN. EQ. N)  G0T0471 
469  CONTINUE 

DO  470  I=NNP1,NP1 
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CAA ( I ) =0 . ODO 
CA ( I ) =0 . ODO 
CONTINUE 

IF(NN.EQ.O)CAAU)=I.ODO 


NOW  THE  FIRST  Hi:  ENTRIES  OF  CA  CONTAIN  THE  S- DOHA  IN  ZEROES  OF  t.'UME 
AND  THE  REMAININGf.HTRIES  ARE  ZEROED  OUT. 


I F (NN . NE. 0 ) CALL  POLCON ( CA , CAA »  0  ,  N ) 


WORK  ON  DENOF.INATOR 


CALL  POLRT(B,TEMP,N,RR,RI,  IER) 

D016  1=1, N 

CR  (I )  =DCMPLX  (RR(I )  ,RI  (I)  ) 

CF ( I ) =1 . 0D00/CR ( I ) 

IF ( I  PR. GE. 2) WRITE (6,1002) 

IFUPR.GE.  2)  CALL  PRCVEC  (CF,N) 
IF(IZTS.EO.O)  GO  TO  900 
D06  1  =  1, N 

CR ( 1 ) = ( -1 . 0 /DELTA) *CDLOG (CR ( I ) ) 
IFUPR.GE.  2)  WRITE  <6, 240) 

FORMAT ( 1  LOGARITHMIC  TRANSFORMATION* ) 
IF (IPR.GE. 1) WRITE (6, 999) 

WRITE (6,2000) 

FORMAT ('  POLES  IN  S  DOMAIN*) 
IFUPR.GE.  1)CALL  PRCVEC(CR,N> 

D03000 1= 1 , N 
CR ( I ) =-CR ( I ) 

CALL  POLCON  (CR,CB,0,N) 


ADJUST  DC  GAIN  CONSTANT 


A2  =  CAA(1) 

B2=CB ( 1 ) 

FAC= (Al/Bl) * (B2/A2) 
DO  603  1=1, NNP1 
CAA ( I ) =CAA ( I ) *FAC 
GO  TO  2010 


DELAYED  PULSE  INVARIANT  TRANSFORMATION 


SHIFTS  NUMERATOR  COEFFICIENTS  FOR  DELAY 


COKT= A ( 1 ) 

DO  300  1  =  1, N 
A (I ) =A ( 1  +  1 ) -CONT*B ( 1+1 ) 

A (NP1 ) =0 . 0 

CALL  POLRT(B, TEMP, N,RR,RI, IER) 

D06 1 1=1 , N 

CR ( I ) =DCMPLX ( HR ( I ) ,RI (I) ) 

CF(I)=1. ODOO/CR ( I ) 

IFUPR.GE.  2)  WRITE  (6, 1002) 

FORMATUX,  'THE  POLES  OF  THE  Z-DOMAIN') 
IFUPR.GE. 2)CALL  PRCVEC(CF,N) 


PARTIAL  FRACTION  EXPANSION 
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D03I=1,N 
C0N1=1 . ODOO 
CON2=0 . ODOO 
D04  J  = 1 ,  N 

cor;  2= cot;  2  *cr  < i )  +a  <n-j+i ) 

IF:i-J) 5,4,5 

5  CO:U=CO!ll*  (l.ODOO-CR(I)  *CF(J)  ) 

4  CONTINUE 

3  CA(I)=CON2/CONl 

c 

c 

C  TRANSFORMATION  OF  DENOMINATOR  AND  NUMERATOR 

c 

c 

224  D02 1= 1 , N 

CR  ( I )  =CDLOG  ( CP.  ( I )  )  /DELTA 

CA (I )=CA (I ) *CR(I) / (I.0D00-CF (I) ) 

2  CONTINUE 

CA(NP1)=0.0D00 

IF ( I  PR. GE. 2) WRITE (6, 241) 

241  FORMAT ( 1  DELAYED  PULSE  TRANSFORMATION’ ) 

I F (I  PR. GE. 0 ) WRITE (6,999) 

226  WRITE (6,1004) 

1004  FORMAT  (  ’  NEGATIVE  OF  THE  POLES  IN  THE  S-DOMAIN  ’  ,  14 ) 
I F  ( I  PR.  GE .  0 )  CALL  PRCVEC(CR,N) 

IFQPR.GE.  2)  WRITE  (6,1003) 

1003  F0PI1AT<  IX,  ’NUMERATOR  CONSTANTS  OF  FACTORIZED  H(S)’) 
IF(IPR.GE.2)  CALL  PRCVEC(CA,N) 

CALL  POL CON (CR,  CB,  0  ,  N) 

D07 1 1=1 , NP1 
71  CAA ( I ) =0 . ODOO 

D09K=1,N 

CALL  POLCON (CR, CF1 , K, N) 

D09J=1,N 

9  CAA(J)=CAA(J)+CF1(J)*CA(K) 

CAA (NP1) =0 . ODOO 
2010  CONTINUE 

DO  450  1=1, NP1 

450  CAA ( I ) =CAA ( I ) +CONT*CB ( I ) 

C 

C 

C 

403  IF ( IPR. GE. 1) WRITE (6,1005) 

1005  FORMAT C  S-DOMAIN  DENOMINATOR’) 

I F ( I  PR. GE. 1 ) CALL  PRCVEC (CB, NP1 ) 

I F ( I  PR. GE. 1 ) WRITE (6,1006) 

1006  FORMAT (’  S-DOMAIN  NUMERATOR’) 

IF  ( IPR.  GE.  1 )  CALL  PRCVEC  (CAA,  HP1 ) 

DO20 1=1 , NP1 
B ( I ) =CB ( I ) 

20  A ( I ) =CAA  ( I ) 

900  RETURN 
END 
C 

END  OF  DATA 
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APPENDIX  B 


A  Tutorial  Example  of  the  Application  of  Eq.  (6) 

Consider  the  setup  of  Fig.  6  where  u  (k)  denotes  the  input  signal  and 
y  (k)  the  output.  The  network  is  known  to  be  of  first  order.  The  measure¬ 
ments  are  made  every  1  ms  for  5  samples,  k  =  0,1,..., 4. 

Suppose  the  following  signals  are  generated  as  a  result  of  a  unit  pulse 
input; 


yO  (k) 

0 

0.5 

0.4 

0.  32 

0.256 

y-\  (k) 

0 

0.  5 

0.9 

1.22 

1.476 

u0(k) 

0 

1 

0 

0 

0 

Uj  (k) 

0 

1 

1 

1 

1 

The  Gram  matrix  of  the 

signals  yQ , 

y-|  and 

°1 

is 

0.57794 

F  - 

1. 37831 

4.7270 

-1.47602 

-4.0960 

4. 

0 

which  yields  the  following  square-roots  of  the  diagonal  cofactors. 

/D  =  1.4597  /D2  =  0.36492  /v  =  0.9123 

Note  that  the  signs  of  these  square-roots  are  chosen  in  direct  correspondence 
with  the  signs  of  the  cofactors  of  the  first  row  of  F  [9].  Now,  substitution 
into  (6)  yields 

(1-0.8  z"1)  Y(z)  =  0. 5  z-1  U(z) 

Clearly,  the  true  parameters  have  been  recovered. 
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y  (k) 
I 


Unknown  network  function 


0-5z 
1  -  0. 8z 


-1 


Fig.  6.  A  first  order  test  system 
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APPENDIX  C 


c 

c 

c 

c 

c 

c  ************************************************************* 

c  * 

C  *  PROGRAM  NAME  '  STOZ  ' 

C  * 

C  *  _ 

c  * 

C  *  S-DOMAIN  TO  Z -DOMAIN  CONVERSION  PROGRAM 

C  * _ 

c  * 

C  *  DEVELOPED  AT  THE 

C  *  UNIVERSITY  OF  SOUTH  FLORIDA 

C  * 

C  * 

C  * 

C  *  _ 

c  * 


C  *  FOR  ASSISTANCE 

C  *  CONTACT 

C  *  DR.  V.  K.  JAIN 

C  *  813/974-2581 

C  ************************************************************* 

C 

c 

S-DOMAIN  TO  Z-DOMAIN  CONVERSION  PROGRAM 

(FROM  A  GIVEN  RATIONAL  TRANSFER  FUNCTION  H  (S)  ,  IT  FINDS 
AN  EQUIVALENT  Z-DOMAIN  TRANSFER  FUNCTION  H  (Z)  BY  ONE 
OF  SEVERAL  METHODS,  E.G.  IMPULSE-INV.  OR  STEP-INV. ) 


GIVEN  THE  CONTINUOUS  DESCRIPTION,  SUBROUTINE  COMPUTES  THE 
EQUIVALENT  DISCRETE  DOMAIN  DESCRIPTION  OF  A  LINEAR 
DYNAMIC  SYSTEM 

STOZ  GENERATES  H(Z)  AND  THE  CORRESPONDING  DIFFERENCE 
EQUATION  FROM  THE  TRANSFER  FUNCTION  H(S> 

THE  INPUT  ARRAYS  A  AND  B  ARE  FILLED  ACCORDING  TO  THE 
DIFFERENTIAL  EQUATION 

B(l)*Y(T)+B(2)*D(l,Y(T))+...+B(N+l)*D(N,Y(T)) 

-A<l)*U(T)-A(2)*D<l,U(T))-...-A(N+l)*D(N,U(T)>  -  0 
WHERE  D (M,  F  (T) )  =  THE  MTH  TIME  DERIVATIVE  OF  FUNCTION,  F 

B  (N+l )  MUST  EQUAL  1 


/ 
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RETURNS  ARRAYS  A  AND  B  CONTAINING  THE  EQUIVALENT  DISCRETE 
DESCRIPTION  STORED  ACCORDING  TO  THE  DIFFERENCE  EQUATION 
B(l)*Y(K)+B(2)«Y(K-l)+...+B(N+l)*Y(K-N) 

-A(l)*U(K)-A(2)*U(K-l>-...-A(N+l)*U(K-N>  =  0 

B  ('. )  ALWAYS  EQUALS  I 

THE  POLES  OF  THE  CONTINUOUS  DOMAIN  MUST  BE  DISTINCT  AND 
NON-ZERO  FOR  THE  TR AN S FORMAT I ON  TO  BE  VALID 


DATA  CARL  SET  PREPARATION 
N  “  ORDER  OF  SYSTEM 

N  (MAXIMUM)  =  ONE  LESS  THAN  THE  DIMENSION  SUBSCRIPT 

ItlTHD  »  0  FOR  THE  IMPULSE  INVARIANT  DESCRIPTION  . 

=  1  FOR  THE  PULSE  INVARIANT  DESCRIPTION  . 

«  2  FOR  THE  TRAPEZOIDAL  INVARIANT  DESCRIPTION  . 

=  3  FOR  THE  LOGRITHMIC  TRANSFORM  DESCRIPTION  . 

IPOLZ  «  1  IF  POLES  AND  ZEROES  ARE  READ 

MUST  BE  COMPLEX  (REAL, IMAGINARY) 

NEGATIVES  OF  POLES  AND  ZEROS  READ;  IE, 

FOR  A  POLE  OF  (S+2) ,  INPUT  +2.0  +0.0 

(2F10.0  PER  POLE',  4  POLES  PER  CARD* 

=  0  IF  DENOMINATOR  AND  NUMERATOR  ARE  READ 
IN  POLYNOMIAL  FORM. 

COEFFICIENTS  ORDERED  FROM  LOW  TO  HIGH  DEGREE, 
DENOMINATOR  INFORMATION  ALWAYS  READ  FIRST. 
HIGHEST  ORDER  DENOMINATOR  COEFF  MUST  BE  1.0 

NOTE;  WHEN  POLE-ZERO  DATA  IS  ENTERED  A  GAIN  CARD 
MUST  FOLLOW  THE  LAST  POLE  CARD. 

IF  THERE  ARE  NO  ZEROS,  USE  BLANK  CARDS  IN 
THE  NORMAL  ZEROS  POSITIONS. 

I ZERO  »  1  ZEROES  OF  Z-DOMAIN  FN  ARE  PRINTED  OUT 


START  EACH  DATA  CARD  SET  WITH  A  DESCRIPTION  CARD, 
CONTAINING  UP  TO  51  CHARACTERS  COLS  2-52 
FIRST  DATA  CARD  CONTAINS: 

N, IMTHD, IPOLZ, IN  3F5  FORMAT,  PLUS  DELTA  IN  F10.0  FORMAT 
SECOND  GROUP  OF  DATA  CARDS  IS  N  POLES  OR  N+l  DENOMINATOR 
COEFFICIENTS.  AFTER  LAST  POLE  CARD,  USE  A  GAIN  CARD. 
LAST  GROUP  OF  DATA  CARDS  IS  N  ZEROS  (OR  BLANKS), 

OR  N+l  NUMERATOR  COEFFICIENTS  (BLANKS  FOR  ZERO  COEFFS) 
THE  DATA  FORMAT  FOR  EACH  OF  THE  SYSTEM  PARAMETER  CARDS 
IS  8F10.0 

AS  MANY  SETS  OF  DATA  CARDS  HAY  BE  RUN  AS  DESIRED 
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c 

C  STOZ  MAIN  PROGRAM 

C 

REAL*8  B ( 20 ) ,A(20) , RR(20) ,RI (20) .DELTA, TEMP (20) 
COMPLEX*16  CR (20) ,CA(20) , CB ( 20) , CAA (20 ) , CAl (20 ) 
1  TEH  ( 20 )  ,  CONI ,  COH 2 ,  CONT 
DIMENSION  TITLE<52) 

C 

C  READ  TITLE  AND  FIRST  DATA  CARD 

C 

100  READ (5,920,END=5999) TITLE 
WRITE (6,910) 

910  FORMAT (6 (/)) 

WRITE (6, 920) TITLE 
920  FORMAT (52A1) 

READ (5,920 ) TITLE 
WRITE (6,930) 

930  FORMAT(/,lX,71( '*’ )  ) 

READ (5, 940 ) N, IMTHD, IPOLZ, DELTA, IZZERO 
940  FORMAT (315, F10. 0,315) 

WRITE (6, 950 )M, IMTHD, IPOLZ, DELTA 
950  FORMAT (//3X, 'N  =', 14 , 5X, 1 IMTHD  = 1 , 14 , 5X, ' IPOLZ 
1  ’DELTA  =' ,G17.10,//) 

NP1=N+1 

NP2=N+2 

NPNP1=N+N+1 

NPNP2=N+N+2 

IF(IPOLZ.NE.l)  GO  TO  300 
C 

C  READ  POLES  AND  ZEROS 

C 

READ (5, 960) (CR (I ) , 1=1 ,N) 

960  FORMAT (8F10.0) 

CALL  POLCON (CR, TEM, 0 ,N) 

DO109  1=1, NP1 
109  B ( I ) =TEM (I ) 

C 

C  READ  GAIN  CARD 

READ (5,960)  RK 
C 

C  READ  ZEROS 

READ (5,960)  (CA (I ) , 1=1 ,N) 

CALL  POLCON (CA, TEM, 0,N) 

DO  209  1=1, NP1 
209  A ( I ) =TEM ( I ) *RK 

GO  TO  310 
C 

C  READ  DENOMINATOR  AND  NUMERATOR  COEFFICIENTS 

C 

300  READ(5, 960)  (B (I) , 1=1 ,NP1) 

READ (5,960)  (A (I ) , 1=1 , NP1 ) 

310  CONTINUE 

C 

C  PRINT  DENOMINATOR  AND  NUMERATOR  COEFFICIENTS 

C 
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WRIT£<6,970) 

970  FORMAT ('  S-DOMAIN  DENOMINATOR') 

CALL  PRVEC<B,NP1) 

WRITE (6,980) 

980  FORMAT ('  S-DOllAIN  NUMERATOR') 

CALL  PRVEC  (A,  tlPl ) 

C 

C  DETERMINE  ORDER  OF  NUMERATOR 
C 

NN-N 

DO  309  I»1,NP1 
II»NP1+1-I 

IF (A ( II ) . NE. 0 . 0 )  GO  TO  400 
309  NN-NN-1 
400  CONTINUE 

WRITE (6, 990)  NN 

990  FORMAT ( 1  ORDER  OF  S-DOUAIH  NUMERATOR  »',I5,//) 
IF(HH.LT.O)  GO  TO  5029 

NNP1-NN+1 

C 

C  FACTOR  DENOMINATOR  TO  FIND  POLES 
C 

IF (IPOLZ. NE. 0)  GO  TO  500 
CALL  POLRT(B,TEMP,N,RR,RI,  IER) 

DO  409  I-1,N 

409  CR(I)«=DCMPLX(RR(I)  ,RI(I)) 

500  WRITE (6 ,991) 

991  FORMAT ('  POLES  OF  THE  S-DOMAIN') 

CALL  PRCVEC  (CR,  N) 

C 

IFdMTHD.NE.  3)  GO  TO  1500 
C 

C  LOGRITHMIC  TRANSFORM 

C 

1000  CONTINUE 

WRITE (6, 9100) 

9100  FORMAT ( / , '  LOGRITHMIC  TRANSFORM') 

C 

C  WORK  ON  NUMERATOR 

C 

IF (NN. EQ. 0 )  GO  TO  1030 

IF ( IPOLZ . NE. 0)  GO  TO  1010 

CALL  POLRT(A,TEM?,NN,RR,RI,  IER) 

DO  1009  I-1,NN 

1009  CA (I ) -DCMPLX (RR(I),RI(I)) 

1010  WRITE <6, 9200) 

9200  FORMAT ( '  ZEROS  IN  S  DOMAIN') 

CALL  PRCVEC (CA,NN) 

DO  1029  I-l.NN 

1029  CA<I)-CDEXP<CA<I) ‘DELTA) 

IF(NN.EQ.N)  GO  TO  1100 

1030  DO  1039  I«NNP1,NP1 
CAA(I)-O.ODO 

1039  CA(I)«0.0D0 

1100  CONTINUE 
C 

C  NOW  THE  FIRST  NN  ENTRIES  OF  CA  CONTAIN  THE 
C  Z-DOMAIN  ZEROS  OF  THE  TRANSFER  FUNCTION,  WHILE  THE 

C  REMAINING  ENTRIES  ARE  ZEROED  OUT. 
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C  WORK  ON  DENOMINATOR 
C 

DO  1129  1=1, U 

1129  CR(I)=CDEXP(CR(I) ‘DELTA) 

C 

C  NOW  CR  CONTAINS  T„J  H  Z-DOMAIN  POLES 

C 

C  FORM  NUMERATOR  AND  DENOMINATOR 

C  Z-DOMAII!  POLYNOMIALS 

C 

IF(NN.EO.O)  CAA ( 1 ) =1 . 0D0 

IF (NM. HE. 0)  CALL  PCSTZ (CA, CAA, 0 , NN) 

CALL  PCSTZ  <CR,CB,0,N) 

C 

C  HOW  CB  CONTAINS  THE  N+l  Z-DOMAIN  DENOMINATOR  COEFFICIENTS, 

C  AND  CAA  CONTAINS  THE  HN+1  NUMERATOR  COEFFICIENTS. 

C 

C  ADJUST  DC  GAIN  CONSTANT 

C 

A1=A(1)/B(1> 

A2=l . 0D0 
DO  1209  1=1, NN 
1209  A2=A2+CAA ( 1+1 ) 

B2=1.0D0 
DO  1219  1=1, N 
1219  B2=B2+CB ( 1+1 ) 

FAC=A1 *B2/A2 
DO  1229  1=1 , NNP1 
1229  CAA ( I ) =CAA ( I ) *FAC 
C 

C  NOW  CAA  CONTAINS  THE  ADJUSTED  Z-DOMAIN  NUMERATOR  COEFFICIENTS 

C  AND  FAC  CONTAINS  THE  GAIN  FACTOR  USED  FOR  THE  ADJUSTMENT 

r~ 

GO  TO  5000 
C 

1500  CONTINUE 
C 

C  NON-LOGRITHMIC  TRANSFORMATIONS 
C 

C  ADJUST  FOR  DIRECT  TRANSMISSION 
C  THIS  ROUTINE  REQUIRES  THAT  B<NP1)  =  1.0 
C 

IF(NN.LT.N)  GO  TO  1510 
CONT=A  (NP1 ) 

DO1509  1=1, N 

1509  A ( I ) =A ( I ) -CONT*B  < I ) 

C 

C  FIND  NUMERATOR  CONSTANTS  FOR  PARTIAL  FRACTION  EXPANSION 
C 

1510  D01529  1=1, N 
CONl=l . 0D00 
CON2=0 . 0D00 
D01519  J=1 ,N 
CON2=CON2*CR ( I ) +A (N-J+l ) 

IF(I-J) 1512,1519,1512 

1512  CONl=CONl* (CR(I)-CR(J) > 

1519  CONTINUE 
1529  CA ( I ) =CON2/CONl 
WRITE (6, 9300) 

9300  FORMAT ( '  NUMERATOR  CONSTANTS  OF  THE  FACTORIZED  H(S)') 

CALL  PRCVEC (CA, N) 
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C  CONVERT  THE  FIRST  ORDER  PARTIAL  FRACTIONS  TO  Z  DOMAIN 
C 

IMTHD-IMTHD+1 

GOTO  (2000,3000,4000),  IMTHD 
C 

C  IMPULSE  INVARIANT 
C 

2000  D02009I-1,N 

CA ( I )  “CA  ( I ) ‘DELTA 
2009  CR ( I ) »CDEXP (CR ( I ) ‘DELTA) 

WRITE (6,3351) 

3351  FORMAT ( 2X, 'Z-DOHAIN  RESIDUES*) 

CALL  PRCVEC  (CA,  N) 

GO  TO  4500 
C 

C  PULSE  INVARIANT 

C 

3000  D03009  1-1, N 

CONl-CDEXP (CR ( I ) ‘DELTA) 

CA (I) «CA (I ) * (CON1-1 . 0D00 ) /CR (I ) 

3009  CR(I)-C0N1 

WRITE (6,3351) 

CALL  PRCVEC  (CA,N) 

GO  TO  4500 
C 

C  TRAPEZOIDAL  INVARIANT 

C 

4000  ICHECK-2 

D04009  I-1,N 
CONl-CDEXP (CR  < I ) ‘DELTA) 

CON2-CA ( I ) / (CR  < I ) *CR  < I ) *DELTA*C0N1 ) 

CONT-CONT+CON2* ( (1 . 0D00-CR (I > ‘DELTA) ‘CON1-1 . 0D00) 

CA  <I)-CON2* (1. ODOO-COU1 ) * (1 . ODOO-CONl ) 

4009  CR( I ) -CONI 
GO  TO  4500 
C 

C  CONSTRUCT  THE  Z  DOMAIN  DENOMINATOR 

C  AND  NUMERATOR  POLYNOMIALS 

C 

4500  CONTINUE 

CALL  PCSTZ (CR, CB, 0 , N) 

DO  4509  I-1,N 
4509  CAA(I)-0.0D00 

DO  4519  K-1,N 
CALL  PCSTZ  (CR,CA1,K,N) 

DO  4519  J«1,N 

4519  CAA  < J) -CAA (J) +CA1 <J) *CA (K) 

CAA (NP1) -0. ODOO 
C 

C  ADJUST  FOR  DIRECT  TRANSMISSION 

C 

DTXC-(0. 0,0.0) 

IF(NN.NE.N)  CONT-DTXC 
CAA ( NP1 ) -CONT*CB ( N  PI ) 

DO  4529  I-1,N 

4529  CAA (I) -CAA (I) +CONT*CB (I ) 

C 

C  SHIFT  NUMERATOR  TO  COMPLETE  PULSE  INVARIANT  TRANSFORM 

C  WHEN  NUMERATOR  HAS  LOWER  ORDER  THAN  DENOMINATOR 
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CONTINUE 

C 

c  - 

c 

c 

C  PRINT  THE  TRANSFORMED  COEFFICIENTS 

C 

5000  CONTINUE 
C 

WRITE (6,9510) 

9510  FORMAT ( 1  POLES  IN  THE  2  DOMAIN') 

CALL  PRCVEC (CR»  N) 

WRITE (6, 9520)  FAC 

9520  FORMAT ( '  GAIN  FACTOR  USED  ='.E14.7,//> 

IF  (IZZEP.O.  EQ.  0)  GO  TO  3323 
NN=N-1 

DO  3321  1=1, N 

3321  B ( I ) =CAA (NP1-I ) 

CALL  POLRT  <B, TEMP, NN , RR, RI , IER ) 

DO  3322  1=1, NN 

3322  CA (I ) =DCMPLX (RR (I) ,RI (I ) ) 

WRITE (6,9530)NN 

9530  FORMAT ( 1  ZEROS  IN  THE  Z  DOMAIN', 14) 

CALL  PRCVEC (CA,NN) 

3323  CONTINUE 
WRITE (6,9540) 

9540  FORMAT  ('  Z-DOtlAIN  DENOMINATOR') 

CALL  PRCVEC  (CB.NPl) 

WRITE (6, 9550) 

9550  FORMAT ('  Z-DOMAIN  NUMERATOR ' ) 

CALL  PRCVEC  (CAA,NP1) 

DO  5019  1=1, HP1 
B  ( I )  =CB  ( I ) 

5019  A(I)=CAA(I) 

1239  FORMAT ( 7  FI 0 . 7 ) 

WRITE (6, 1239) (B(l) ,I=l,NPl) 

WRITE  (6,123  9)  (A  (I )  ,  1=1  ,NP1 ) 

GO  TO  100 
C 

5029  WRITE  <6 , 9560) 

9560  FORMAT  (/,  IX,  'NUMERATOR  ORDER  LESS  THAN  ZERO',//) 
5999  STOP 
END 
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APPENDIX  D 


GEORGE'S  THEOREM 


This  theorem  helps  evaluation  of  the  quadratic  volterra  response  (or 
the  bilinear  response)  by  inspection  from  the  associated  response. 

Theorem  If 

1 


^ ( S1 *So)  “ 


C(s,+  s_) 


(2) v  1*  2'  (Sl+a)(s2+b)  1  2 


(Bl) 


then 


Y2(s) 


C(s) 


s  +  a  +  b 
Proof:  It  is  readily  shown  that 

-jco 


*2(s) 


1 

2*j 


-J 


C(s) 


(s^+a) (s-s^+b) 


ds^ 


(B2) 


(B3) 


(see  Bush  [11]).  Then  the  desired  result  follows  immediately  by  use  of  the 


residue  theorem. 
Corollary  Let 


Y(2)(S1,S2) 


-as^  -Bs2 

e _ _e _ 

(sL+a) (s2+b) 


C(s) 


(BA) 


If  8  >  a,  then 


Y2(s) 


-(B-a)a  e~BsC(s) 

(s  +  a  +  b) 


(B5) 


If  a  >  B  then 


-b(a-B)  e  asC(s) 

(s  +  a  +  b) 


Y2(s)  =  e 


(B6) 
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